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ABSTRACT 


Explicit difference equations are presented for the solution of 
a signal of arbitrary waveform propagating in an ohmic dielectric, 
a cold plasma, a Debye model dielectric, and a Lorentz model dielectric. 
These difference equations are derived from the governing time- 
dependent integro-differential equations for the electric fields by 
a finite difference method. A special difference equation is 
derived for the grid point at the boundary of two different media. 
Employing this difference equation, transient signal propagation in 
an inhomogeneous media can be solved provided that the medium is 
approximated in a step-wise fashion. The solutions are generated 
simply by marching on in time. By appropriate choice of the time and 
space intervals, numerical stability and convergence are always ob- 
tained. Numerous examples are given to demonstrate the wide range of 
applicability of the difference solution. These include: the trans- 
mission and reflection of an electromagnetic pulse normally incident 
on a multilayered ohmic dielectric; a step-modulated sine wave 
propagating in a dispersive media, a problem originally considered 
by Sommerfeld and Brillouin; the reflection of a short gaussian pulse 
normally incident on an inhomogeneous lossy cold plasma with a 
longitudinal d.c. magnetic field, and many others. It is concluded 
that while the classical transform methods will remain useful in 
certain cases, with the development of the finite difference methods 
described in this dissertation, an extensive class of problems of 
transient signal propagating in stratified dispersive media can be 
effectively solved by numerical methods. 
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CHAPTER I 
INTRODUCTION 


The propagation of a transient signal in dispersive media is an 
important and interesting problem in electrical science. Take, for 
example, the historically famous problem considered by Sommerfeld[l]. 

In this work he showed that the front of a signal propagates exactly 
with the speed of light in free space regardless of what the group 
velocity might be. Thus it supported, at a controversial time[2], 
one of the fundamental postulates in Einstein's special theory of 
relativity. More recently, there has been an interest in the remote 
sensing of the thickness of ice layers over water, and in the sensing 
of subsurface geological structure using short pulse techniques[3]. 

Other applications are the distortion of the waveform of a signal 
reflected from the ionosphere[4], and the measurement of the electrical 
properties of dielectrics utilizing time-domain reflectometry[5] . 

To analyze these transient problems, Fourier and Laplace transform 
methods are usually employed. Ahtough these standard methods apply to 
many cases, their failure to accommodate others has warranted additional 
and new methods of solution. In the case of cold plasmas, for instance, 
Bowhill[6] proposed the multiple-scattering technique, and Field[7] 
suggested the method of characteristics. This dissertation presents 
some finite difference methods for the solution of transient signal 
propagation problems in stratified dispersive media. 

The numerical solution of partial differential equations has its 
root deep in the past. In 1928, Courant, Friedrichs and Lewy[8] pub- 
lished a celebrated paper in which they proposed originally the idea 
of solving the wave equation, the diffusion equation and the Laplace 
eqi ation by algebraic operations. In particular, when they replaced 
thvi wave equation 


( 1 ) 


2 2 

^ - c^ ^ = 0 
2 c 2 

at az 


by the difference equation 


E(z,t+At)-2E(z,t)+E(z,t-At) _ 2 E(z+Az ,t)-2E(z ,t)+E(z-Az ,t) _ ^ 

p . C p U j 

at^ az^ 


they found that the numerical solution always converges if the grid 
ratio cAt/Az is less than or equal to' one. Furthermore, if the 
equality sign is chosen, the resulting difference equation 


1 


■4 
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(2) E(z,t+At) = - E(z,t-At) + E(z+Az,t) + E(z-AZ,t) 

generates exactly the analytic solution[9]. More than 30 years later, 
in 1970 Chiu[10] carried this remarkable solution through an inter- 
face boundary for the first time. He derived an appropriate difference 
equation for this case and showed that the numerical method again 
produces an exact solution. 

We shall extend these techniques to the integro-differential 
equation of the type: 


(3) 



c^ ^ + a ||-+ bE(z.t) 


K(t-g)E(z,6)dg. 

^0 


As will be seen this type of equation, or various special cases of 
it, arises in a wide variety of problems when a plane electromagnetic 
wave propagates in linear dispersive media. 

In Chapter II we shall discuss in detail the numerical solution 
of (3) in a homogeneous medium. A procedure will then be presented 
to derive the appropriate difference equation for a stratified medium. 

As an example the wave equation is considered here. By approximating 
the inhomogeneous dielectric medium with a large number of thin layers, 
it will be shown how one of the most challenging boundary-value 
problems in electromagnetic field theory, the reflection and trans- 
mission of a plane electromagnetic wave from an inhomogeneous di- 
ele:tric slab[ll], can be solved as an initial-value problem - by 
Uv-ilizing the step-modulated sinusoidal wave as the incident signal 
and marching the time domain solution long enough to allow the 
steady-state to be established. 

In Chapter III a difference equation for a stratified lossy 
dielectric medium is first derived. This difference equation is then 
utilized to obtain the waveform of a unit-step signal propagating in 
a homogeneous lossy dielectric medium, the transmission of a gaussian 
pulse from the air into the earth, and the reflected waveform of a 
sine-squared pulse from a three-layered lossy dielectric. In all 
cases normal incidence is assumed. 

In Chapter IV a numerical solution to Sommerfeld's "About the 
propagation of light in dispersive media"[l] is investigated, A 
difference equation is derived from the governing integro-differential 
equation for the electric field. Employing the resulting difference 
equation, we have been able to observe the real-time propagation process 
vividly on a high-speed digital computer having a CRT display facility. 
The dynamic behavior of the precursor of a signal thus can be seen. 

We shall reexamine this historical problem and discuss some of our 
findings. 



Chapter V deals with transient signal propagation in a stratified 
cold plasma in which a constant magnetic field parallel to the direction 
of propagation can be assumed. The coupled integro-differential 
equations for the electric fields are first derived from Maxwell's 
equations and the equation of motion for the electrons. They are next 
transformed into an explicit difference equation; the solution is 
then obtained simply by marching on in time. By approximating the 
plasma medium with a large number of thin layers, it is shown that the 
time history of a signal of arbitrary waveform in an inhomogeneous, 
anisotropic {due to the presence of a constant magnetic field), lossy, 
cold plasma can be easily obtained. Again, only the case of normal 
incidence is considered. 

In Chapter VI we calculate the reflected waveforms of a unit- 
step signal normally incident from air to water directly as well as 
with an ice layer on it. Here the water is considered as a Debye 
dielectric[12]. We first derive the governing time-dependent 
integro-differential equation for the electric field propagating in 
this type of dispersive medium, and then deduce the corresponding 
difference equation. The solution is then generated again by 
marching on in time. 

Finally in Chapter VII we summarize the results obtained, and 
discuss briefly the application to nonlinear dispersive media as de- 
veloped elsewhere[13,14], a subject of great interest in recent years[15]. 

It is concluded that while the classical Fourier and Laplace 
transform methods will remain useful in certain cases, with the de- 
velopment of the finite difference methods described here as well as the 
a/ailability of a high-speed digital computer, an extensive class of 
problems of transient signal propagation in stratified dispersive media 
can be effectively solved by numerical methods. 



CHAPTER II 

DIFFERENCE METHODS FOR AN INTEGRO-DIFFERENTIAL EQUATION 


Although the numerical solution of partial differential equations 
by finite difference methods is well-established[16,17,18], little is 
known of its application to the initial -value problem of the telegraph 
and related equations[19]. In this chapter we shall discuss the 
numerical methods for the integro-differential equation (3) in detail. 
We shall start from the wave equation, then add one lower order term 
after another. In this manner, the lossy wave equation, the telegraph 
equation, and the integro-differential equation itself are treated. 

We shall also present a procedure to derive the difference equation 
for a grid point at the boundary of two different media, in order to 
be able to treat propagation in inhomogeneous media. 


A. The Wave Equation 

As was mentioned earlier, if we replace the wave equation 

9t az 


by the difference equation 


E(z,t+At)-2E(z,t)+E(z.t-At) _ 2 E(z+Az,t)-2E(z,t)+E(z-Az,t) 


At‘ 


Az 


or 


(4) 


E(z,t+At) = -E(z,t-At) +2 

\2 


-f 


cAt r 


E(z,t) + 


AZ j 

+ (^] [E(z+AZ,t) + E(z-AZ,t)] 
then the numerical solution always converges if the stability criterion 


(5) 


7 ^ <1 
AZ - 


is satisfied. Furthermore, if the equality sign is chosen, i.e. 


4 



5 


(6) AZ = CAt, 


then the resulting difference equation 


(7) E(z,t+At) = - E(z,t-At) + E(z+Az,t) + E(z-Az,t) 


generates exactly the same values as the analytic solution. This is 
to say that as far as the numerical values of the solution are con- 
cerned, one could not distinguish whether they were either obtained 
from solving the intial value problem of the wave equation or ob- 
tained from (7) by marching on in time. This fact has long been 
known by numerical analysts[9,17]. 

In order to visualize this remarkable property, let us consider 
the Taylor series expansions 

9^E _ E(z+Az,t)-2E(z.t )+E(z-Az.t) 3^E az^ 

.,2 " ._2 "“4 1 ^ **• * 


A similar expression can be obtained for the time variable t. Now 
if we approximate the second partial derivatives by only the first 
term on the right-hand side, obviously a truncation error results. 
Therefore the replacement of the wave equation by (4), which is now 
written in the notation E^ = E(iAz,nAt) as 

( 1) Ef ’ = -E?-'' + 2 1 

introduce a truncation error 



(9) 


T 


1 

12 



1 dh 

C^ 9t^ 


(cAt)' 


higher order terms. 


Equation (8) tells us that if the initial data at the time steps 
n= 0,1 are given, then the fields at subsequent times can be ob- 
tained recursively. However the result will be correct only if a 
proper value of cAt/Az is used. 


A procedure to determine the permitted value of this ratio, or 
the stability criterion as it is called in numerical analysis, is 
described by Hildebrand[17, pp. 235]. Let the solution of (8) be 
the form 





( 10 ) 
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where a,e are constants, with a real and j =\/TT. Since the number 
of time steps is arbitrary, the index n can be increased without 
bound. Therefore the magnitude of 6 in (10) must be less than one 
for a finite solution. This condition is now used to determine the 
stability criterion. 


Substituting (10) into (8) yields the secular equation 


( 11 ) 


- 2 


1 - 2 


( CAt 


1 AZ 

) 2 


a + 1 = 0. 


I 6- = 2 

i=l ^ 


’ - 2 n? 2 


where 3i , i=l ,2 are the roots of (11), the requirement of la. I < 1, 
i=l ,2 implies ^ 


1 




• 


This inequality is true if and 
bracket is less than one. 


cAt 

AZ 


F 


• 2 a 
Sin j £ 


1 . 


only if the absolute value inside the 


Since a is real, the value inside the absolute sign is always 
positive, 


This inequality is true if and only if 

( 12 ) 1 

for all real a.^ This is the well-known stability criterion for the 
numerical solution of the wave equation by the finite difference 
method[16]. 

If the equality sign in (12) is chosen, three remarkable properties 
are seen: 
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1. The time increment At is the largest one permitted, 
thus allowing a problem to be solved with least number 
of time steps. 

2. Calculations at each step are reduced to a minimum because 
(7) is the simplest version of (4). 

3. But above all, exact solution is obtained from the numerical 
method because the truncation error in (9) has now vanished. 

It is indeed remarkable all these sought-after properties in a 
numerical solution of differential equations happen at the same time. 

We have therefore identified the scheme Az = cAt optimum in every 
sense for the numerical solution of the wave equation by the finite 
difference method. This scheme seems to force the difference equation 
to march along the characteristics of the wave equation. Since an 
initial disturbance propagates along the characteristics of the wave 
equation without changing form, it is not surprising that the dif- 
ference equation produces an exact solution. 

The above conjecture can be verified by showing the difference 
equation satisfying the general solution of the wave equation. To 
do so, let the general solution of the wave equation be the form 

E{z,t) = f^(z-ct) + f2(z+ct). 


We compute and verify that 


-E(z,t-At) + E(z+Az,t) + E(z-Az,t) 

= -[f-|(z-c(t-At))+ f 2 (z+c(t-At))] + 

+ [f-| (z+Az-ct)+f2(z+Az+ct)] + [f-| {z-Az-ct)+f2(z-Az+ct)] 
= f.| (z-Ct-Az)+ f^Cz+Ct+Az) 

= f^ [z-c(t+At)]+ f 2 [z+c(t+At)] 

= E(z,t+At) 

Q.E.D. 

In these manipulations the relation Az=cAt has been used. 

We now proceed to extend these results to a number of more 
general situations. 
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B. The Lossy Wave Equation 

The numerical solution to the lossy wave equation 

(13) 

n 8z 


is also known [17, pp. 311]. While the first two terms are approxi 
mated, as usual, by (8), the first derivative is approximated by 
the central finite difference 


Eu,t+At)-h('z.t-Atl ^ "i 
^ 2At ^ ~ 


cn+l ^ rn-1 


for a better accuracy. Equation (13) becomes 


-fTisr- 

2 1 , 

+ fcAtf (^n + gn > 

\azJ ^S‘+1 i-Kj . 

Its stability criterion can be determined by following the same steps 
as (10)-(12) to obtain 

2 ''■2 f^fsin^ 2. 1 - 

(16) _ 2 ^ + ^=0 

1 + ^ 1 + ^ 

The requirement for [ 3 i| £^l,i=l,2 leads to the stability criterion 
for the lossy wave equation 

,17) £|t</T71f _ 

Since we shall consider only a > 0, (17) is always satisfied if the 
scheme Az=cAt is selected. Consequently (15) is simplified to 


a^t ] rn-l . pn pn 
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which is an appropriate difference equation for the numerical solution 
to the lossy wave equation by the finite difference method. 

C. The Klein-Gordon Equation 

The addition of another kind of lower order term to the wave 
equation yields the Klein-Gordon equation[20] 

(19) -M- - c^ + bE = 0. 


Although a numerical solution to this equation has been discussed[19] , 
it has never been carried out according to the best knowledge of the 
author. Nevertheless, at the first glance the difference solution 
for this case seems quite straight-forward. One simply adds the lower 
order term to (8) and obtains 


(20) = - e""^ + 2 






(e" +e" 
^^i+1 S'-l 


)-bAAj. 


However, the stability criterion for this case is found to be 


(21) Jl -ibst^. 


Hence the scheme Az=cAt can no longer be chosen if b>0. This is the 
case for cold plasmas but not for the Debye dielectrics as will be 
seen. 


One may of course select a scheme such that (21) is satisfied 
and proceed to obtain a numerical solution. In so doing not only 
have all the remarkable properties of the difference solution for the 
wave equation been lost, but as will be seen later, the scheme AZ=cAt 
is found to play an essential part in the practical aspects of the 
numerical solution. There is another reason to reject the above 
idea however. It was conjectured previously that the scheme 
Az=cAt is not just a numerical convenience, but it also coincides 
with the characteristic theory of the hyperbolic partial differential 
equation, in that this scheme forces the difference equation to march 
on the characteristics of the wave equation. Since the presence of a 
lower order term in the wave equation does not affect its character- 
istics[21], it is anticipated that the optimum scheme for the wave 
equation should hold true for the Klein-Gordon equation. It turns out 
this is correct. 
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For the time being, let us accept the fact that if we approxi- 
mate the electric field in the lower order term by three con- 
secutive times such as 


(22) e; . E" + i E'J"' 4 e;-i] _ 

then it can be shown that the stability condition for the resulting 
difference equation 


(23) 



1 

1 + \ bAt^ 



fcAt 



)J 


+ 
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agrees exactly with that of the wave equation. Consequently the optimum 
scheme Az=cAt can be selected againl Equation (23) then simplifies to 


(24) 


pn+1 _ r-n-1 

ti - - h. 


+ 


1 + 1 bst^ 







w. ich is an appropriate difference equation for the Klein-Gordon 
equation. We did not find out the scheme (22) by chance however. In 
the next section we shall demonstrate how it was derived. 

D. The Telegraph Equation 

Combining the lower order terms in the two previous cases, we have 
the telegraph equation 


(25) LE _ ^2 ^ ^ ^ ^ ^ ^ 

When a=GL + RC, b=RG, and c^ = 1/LC, (25) is the familiar transmission 
line equation in which R, G, L, C are the resistance, conductance, 
inductance, and capacitance per unit length of the transmission line. 
Our goal here is to show that the electric field in the last term 
should be approximated by three consecutive times as suggested in 
(22). The approach will be to find a numerical solution which is 
consistent for both the second order partial differential equation 
and the first order equations from which it was derived. 



n 


To deduce this result, let us return to the second order wave 
equation 

— ' c — 2 - 0. 

8z 

It has been shown that if the scheme 
(26) Az * cAt 

is selected, the resulting difference equation 


(27) 


Ef = - . E",^ . e;.^ 


generates exactly analytic solution. This must be true if the cor- 
responding first order system of equations, the Maxwell's equations 


3H 


ar * sf" 


(28) 


^0 


9t 


8E 


9Z 


= 0 , 


are used as the starting point for the same problem. Therefore, if we 
enp.'oximate (28) by the difference equations[16 , pp. 262] 


(29) 




un+1/2 mTH-1/2 

”i-M/2 ' 1-1/2 

Az 


= 0, 


nn+3/2 i,n+l/2 

^i+1/2 ' ^'-H/2 


+ 


s'?!] - 

1+1 T 

Az 


= 0 


or 


(30) 


Eti'*’! = E*) _ 

i i 


At fL.n+1/2 nn+l/2x 
G Az ^1+1/2 ■ "i-1/2^’ 


mO+ 3/2 ^ nn+1/2 At /pn+1 
^+1/2 ^+1/2 ' y^Az ^1+1 



). 
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we should be able to obtain (27) by eliminating the magnetic field 
in (30). Indeed, we substitute the second difference equation in 
(30) to the first difference equation and obtain 


E^+1 = E? - ^ 
1 1 


At 


= E" 

^i e_Az 


_ At ,pn _pnJ Ln 
Vl/2 u^Az ^^i+1 ^i^ fi 

/|jn”l/2 ijn”l/2> , /pn j. rn \ 

'"l+l/2 * '^i+1 ■ ^^i+1 ^ 


l/2_ At /pH pn xi) 

1/2 V? 


= eJ + (eJ-eV-'') + (E^^, - 2E" + E^_,) 


= - Ef ’ - E^,i + e;_, 


which is exactly identical to (27). 

Now we can do the same thing to the lossy wave equation 


(31) 


J ^+o_ ii _ n 
..2 .2 ^ 3 t “ 

3t 92 0 


The difference equation for this lossy wave equation has been given 
in (18). It is 


(32) 


^n+1 


1 


1 + 


gAt 

2e. 


y Ze^j^i ^ ^-+1 ^ ^i-1 


The corresponding Maxwell's equations 


(33) 


3Ex 9H 

e vr~ + + gE =0, 

0 9t 9z X 


9H„ 9E 

u — ^ + — — = 0 

^0 3t 9z 


can be written in difference equation form 
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(34) 


or 


(35) 


-n+1 


- E'. 


un+1/2 _ un+1/2 

f 1 / o n * 


At 


'i-H/2 


Az 


+ e") = 0, 


nn+3/2 nn+l/2 nll + l pH+1 

„ "i+1/2 ■ ”i+1/2 4 . i+1 ■ i 
*^0 At Az 



1 

1 + 



cjAt A pH At 

2e^ I 1 "e^Az 

0 / 0 


/□n+1/ 2 
^^i+1/2 


un+1/2 
■ ^1-1/2 




„n+3/2 _ un+1/2 
"i+1/2 ”i+l/2 


At / pn+1 rn+1 X 

y Az ^1+1 " 1 > 

0 


Notice that in (34), the electric field in the lossy term has been 
approximated by two consecutive times. It is desired to show that 
(32) can be deduced from (35) by eliminating the magnetic field. 


We first compute 

(36) 

At /un+1/2 un+1/2 

s^Az ^1+1/2 " 1-1/2) 


At un-1/2 _ At /pn _pnx 
£ Az ”i+l/2 vt Az '^1+1 ^i ^ 

O L 0 



n-1/2 

i-1/2 


At 


y.AZ ^i-1^ 


At /un-1/2 

£ Az ^1+1/2 
-0 


un~l/2\ / pO opn.pn ^ 

^•-1/2^ ^S*+r^^i^^i-i^ 





gAt j-n-1 

2^0 j 


(e" 

^1+1 


- 2E" + E*? , ) 
1 1 -r 


and substitute it to the last bracket in the first difference equation 
in p5) to obtain (32). The important feature in (34) is the approxi- 
mation made in the electric field in the lossy term. Thus, for the 
transmission line equations 


C 



9i . P 


0 , 




Ri 


= 0 . 


(37) 



we can approximate the voltage v in Gv and the current i in Ri by 
its average values at two consecutive times as 
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n+l_ n -n+1/2 _ -n+1/2 


At 


•n+3/2 .n+1/2 

I j I T y o I • 


n+1 _ n+1 

L .'j +V^ 'j+1/2 j+1 + p i /.n+3/2 . .n+l/2v 

At ^ Az ^ 2 ^^j+1/2 ^ 'j+1/2^ 


* 0 


or 


(38) 


+ [l - 5^]v') - (j"+l/2,. n+1/2, . „ 

2C f: 1' 2C p CiZ ' J+1/2 'j-1/2' 


fi + Mt|. n+3/2 RAtV-n+1/2 At /,n 

2L /’j+1/2 • 2rj'j+l/2 ■ UF <''j 


J+1 j ' 


Now we eliminate the current i in (38) to obtain a difference equation 
for the voltage v by computing first the term 



j 


and substituting it into the first difference equation in (38) to obtain 
a difference equation for the voltage 


(39) 


1 + 


(GL+RC)At RG(At)^ n+1 
2LC 4LC ''j 


1 


(GL+RC)At 

2LC 


, RG(At)^ 
4U: 



+ 




RG(At)^ 

4LC 


v 


n 

j 


in which the scheme 
(40) Az At 


has been used. 
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It can be verified easily that if we replace the second order 
telegraph equation 


(41) 


LC 


^ ^ + (GL+RC) Iy +RG V = 0 


3t 3z 
by the difference equation 

(42) LC - 

At 






,, -2v^+v'? , 

J— + (GL+RC) J ^ 

Az'^ 


and use the scheme (40), then (39) will result. 


It is seen that the voltage v in the last term RGv in (41) must 
be approximated by its value at three consecutive times in order that 
the difference equations for the first order (transmission line) and 
the second order (telegraph) partial differential equations are 
compatible. It was from here that we have identified the scheme for 
the Klein-Gordon equation as presented in (22). 


E. The Source Term 


We come finally to consider the effect of adding an integral 

S(z,t) = f K(t-e)E(z,6)d3 
"^0 

to the wave equation to produce the integro-differential equation of 
Eq. (3). It is always possible to reduce the integral to recursive 
form by splitting it into two parts 

ft-At ft 

S(z,t) = K(t-3)E(z,3)dS + K(t-3)E(z,3)d3 

”'o -^t-At 

and integrating the last integral from t-At to t by the trapezoidal 
rule of integration. The result is given by 

ft-rAt 

(43) S(z,t) = j K[(t-At)-3+At]E(z,3)d3 + 

Jo 

+ [K(0)E(z,t)+K(At)E(z,t-At)]. 
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Here the first term can always be written in terms of the functional 
S(z ,t-At} . 

So far, we have examined the numerical solution of (1) in a 
homogeneous medium by finite difference methods in detail. By 
recognizing the optimum scheme Az=cAt for the wave equation, we have 
concluded that this scheme can be chosen for the lossy wave equation, 
the Klein-Gordon equation and the Telegraph equation. The corresponding 
difference equations for each of these cases have been gi ven explicitly 
in (18K (24), and (39). The stability for each of these difference 
equations has also been analyzed. On the other hand, we have not been 
able to analyze the stability of a difference equation resulting from 
an integro-differential equation. Stability analysis in this case is 
not available at the present time. Under such circumstances, one has 
recourse to intuition and numerical evidence. Fortunately, in all 
our problems considered later, we have never experienced numerical 
instability of any kind. 

The developments up to this point are strictly valid for homogeneous 
media. Their extensions to stratified media are now in order. In the 
next section we shall consider the wave equation and show how the 
propagation of an electromagnetic wave in an inhomogeneous dielectric 
medium can be solved by the difference method. The other applications 
will be followed up in subsequent chapters where some interesting 
transient problems in dispersive media will be discussed. 

F. Difference Equation for a Stratified Dielectric Medium 


To extend previous techniques to a stratified medium, it is 
worthwhile to recall again the remarkable property of the scheme 

(44) Az = CAt 

for the numerical solution of the wave equation by the finite difference 
method, which is that the resulting difference equation 


(45) 


= - e'?-"' 


+ 

“"i+l 


+ E 


n 

i-1 


generates an exact solution. Although this fact has been known for a 
long time[9,17] it was believed to be of little value for practical 
purposes because the analytic solution for the initial -value problem 
of the wave equation is very simple[22]. In 1970 Chiu[10] in the study 
of stress wave propagation in an elastic bar with discontinuities, made 
an important contribution by carrying this remarkable difference solution 
through an interface boundary for the first time. We shall see this can 
also be done for a multitude of interface boundaries. 
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The generalization of the scheme Az=cAt to a problem involving 
an interface boundary is not difficult to visualize. Let us consider 
this situation as shown in Fig. 1. If we subdivided each region such 
that 


(46) 


Az^ = C^At, 
Az 2 = C2At , 






the difference equation (45) then applies to all interior grid points 
except the interface grid point. In order to complete the solution, 
an appropriate difference equation for this interface grid point must 
be found. This can be accomplished by considering the following 
mathematical procedure. 




2" 

A i 

z* 

^ 

w ^ 



Fig. 1— Interface between two dielectric media. 


Referring to Fig. 1» we see that the wave equation is 
defined everywhere within each region, hence it is valid in the 
vicinity of the interface boundary. Therefore, the wave equation at 
and z" respectively are 


2 a^E(z~^,t) _ 0 

..2 ^2 .2 ^ 

9t 9z 

_ 2 9.^E(z~.t) _ 

.*2 ""l ,_2 ^ 


( 47 ) 
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At the interface boundary, the boundary conditions on the continuity 
of the tangential electric and magnetic fields, lead to the con- 
dition 

E(z"‘’,t) = E(z",t), 

(48) 

3E(2~^.t) _ 3E(z~,t) 

9z 3Z » 


the second of these being a consequence of the fact that the time 
derivative of the magnetic field is also continuous. 


To derive a difference equation from the mathematical problem 
as posed in (47) and (48), we expand the fields at z'^’+AZp and 
z'-Az-j respectively by the Taylor series 

E(z--Az,.t) = E(z-.t) - AZ, + "A 

(49) ^ o . _2 


E(z‘^+Az^,t) = E(z‘^,t) + AZc> + 




9Z 


in which the higher order terms have been neglected. The second partial 
derivatives with respect to z in (49) are first substituted by (47) to 
obtain 


(50) 


E(z“-Az^ ,t) 


= E(z-,t) - AZ, + 


AZ. 


ar 2c 


2 * 


E(z%Az,,t) = E(z‘*',t) +^^^f^AZ2 + ^ 

' ar 2c^ 


We eliminate the first derivatives by the use of the boundary conditions 
in (48) and then replace the second partial derivatives with respect 
to time t by the central finite differences 


3^E(z-,t) _ E(z-,t+At)-2E(z-.t)+E(z-,t-At) 

3t^ At^ 

Doing all these steps, we finally obtain a difference equation for the 
interface grid point 
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(51) E(z ,t+At)= - E(z ,t-At)+ — C/TjE(z -Az^ ,t)+yF 2 E(z^+Az 2 ,t)]. 

in which the condition (46) has been employed. 

Since (51) is derived from a canonical problem, it applies to 
any interface boundary and hence a multiplicity of plane boundaries 
as well. Furthermore, it reduces to (45) when the dielectric medium 
at both sides of an interface grid point i is the same. Consequently, 
(51) can be used to obtain the time history of a signal in an in- 
homogeneous dielectric medium provided that the medium is approxi- 
mated in a step-wise fashion. 

To be explicit, the general difference equation for a multilayered 
lossless dielectric medium is given by 


(52) 




+ 


^?-l 



where ei ,g 2 are the permittivities of the medium immediately to the 
left and to the right of the grid point i. Note that (52) must be 
used in conjunction with the modeling 


(53) Az^ = ■ - - At 

foi the thin layers. The subscript i denotes the ith thin layer and 
At is arbitrary. Because of the constraint in (53), exact modeling 
for the slab is possible only when At is extremely small. In practice, 
high precision is seldom required, and the approximation (53) produces 
results of adequate accuracy. 

G. A Numerical Example -- Reflection of a Step-Modulated Sine 
Wave from an Inhomogeneous Slab ~ 

We illustrate the validity of (52) by the example shown in 
Fig. 2. Here a step-modulated sine wave is normally incident on an 
inhomogeneous dielectric slab with the profile 


(54) e(z) = i (e^+Cj) + 1 (c^-Ejjcos f- 

The slab is first divided into thin layers according to (53), (52) is 
then used to produce the reflected waveform of the incident signal. 
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Fig. 2— Reflected CW wave from an inhomogeneous 

dielectric slab versus N, t=N^t , a)At=2ir/36. 


L= 


2k d 

a 


1,2 „ 1 2 


(ei+e^), ^ - 2. 
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Figure 2c shows the real-time reflected waveform as observed 
at the air-dielectric interface boundary for the case of a normalized 
slab thickness L=0.75 versus time step N. The steady-state value of 
the reflection coefficient obtained by Casey[ll] is shown in Fig. 2b. 
Vje see that the amplitude of the time domain solution at late time 
does agree with the frequency domain solution. According to Casey's 
results, there shall be no reflection of electromagnetic energy at 
steady-state when L=1.650 and 2.825. Indeed, Fig. 3 confirms this 
result. By varying the frequency of the incident signal and repeating 
the calculations for each case, the magnitude of the steady-state 
reflection coefficienct can be obtained. By comparing the steady- 
state waveform of the time domain solution with a reference signal, 
phase information can also be available. 

It is interesting to see that the difficult boundary-value 
problem of plane electromagnetic wave propagation in an inhomogeneous 
dielectric medium[23,24] can be solved in this fashion. In a single 
calculation, the reflected field, the transmitted field, and even 
the fields in the medium are all available during both the transient 
period and at steady-state. However, only the case of normal in- 
cidence is considered here. 

H, Some Practical Aspects of a Difference Solution 

In this section we shall discuss various aspects of the numerical 
solution. They include the simulation of the infinite boundary, the 
moving time window calculation technique, proper treatment of initial 
conditions, numerical modelling, etc. 

To illustrate these techniques we consider the problem of a 
jhort pulse incident on an inhomogeneous dielectric slab of finite 
thickness (see Fig. 4). This slab is divided into thin layers such 
that 



where Az^.e^; i=2, 3, 4, ••• are determined from e{z). 

A technique to accomplish this is to determine a temporary thick- 
ness 


AZ 


2 


At 


and obtain the permittivity of layer 2 by 
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(56) ^2 " f • 

The actual thickness of layer 2 is then determined 


AZ 


2 



The above process is then repeated (considering z=az 2 as the initial 
guess for Az' 3 ) to determine 03 , AZ 3 ; e^, Az^ ••• and so on. In 
general, it is impossible to model the slab exactly unless t is 
extremely small. However such high precision is usually not re- 
quired and the problem can always be solved approximately. Once 
the permittivity in each layer is determined, (52) can be used to 
generate the fields at all grid points. 

Suppose the generator at z < 0 (see Fig. 4) is turned on at n= 0 
As the time step marches on, the signal arrives at the air/dielectric 
boundary; reflection and transmission occur. The reflectes wave 
will eventually reach the send-end boundary at z<0. Physically this 
reflected signal could propagate through the boundary to the left and 
never enter the picture again. However, this is not the case 
numerically. The zero value of the short pulse at later times 
acts like a perfectly conducting boundary. So the reflected signal 
from the slab will appear to be re-reflected from the sending-end 
back toward the slab. Here we see the numerical solution generating 
a non-physical signal. This fact applies to the far-end boundary as 
well. Unless some kind of numerical technique is incorporated at 
both boundaries, a large number of grid points is needed to model the 
problem so as to prevent these non-physical signals from coming back 
to the observation point. Consequently the usefulness of the dif- 
ference method will be severely limited. 

For the far-end boundary, it turns out that these problems can 
be avoided in a simple way by terminating the infinite boundary at 
a point just outside the slab. The value at this grid point, rather 
than being computed from the difference equation, is assigned by 


fC’7\ z: 

*^IEND ‘^IEND-1 . 

This follows from the fact that a signal at the point Eiemd_| 
propagate without distortion to if* space. Hence (57) 

simulates the infinite boundary exactly in this case. A similar 
technique exists for sending-end boundary. If we put the generator 
at i=l and the air/dielectric interface at i= 2 , then when we compute 
the fields at d (see Fig, 4) from (52) 


2 



<^0 ^"-1 



(58) d = -b + 
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ig, 4- -A signal normally incident on an inhomogeneous 
dielectric slab in space-time grid points. 
z=iAz, t=nAt. 
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It is important to interprete E^_-j correctly. Physically, the value 
at each grid point represents the total field at that point. Therefore, 
the quantity E!?_-| has two components: an incident field to d gen- 

erated at c, and a reflected field from b. Thus, 

E'?_i = incident field + reflected field from b, 

(59) 

= c + (b-a) 


Once the value of d is determined, the reflected wave can be picked 
up at i=2 by 

(60) E*^(n t) = d-c 

In this way all the unnecessary calculation in the free space has been 
eliminated and generation of a non-physical signal has also been 
avoided. It is worthwhile to mention that all these techniques are 
possible because the scheme Az=cAt has been selected in the difference 
solution. 


When the far-end boundary is other than free space, one can always 
terminate the infinite boundary at the point 


(61) I 


END 


= I 


<"->cbs> 


obs 


+ 1 


(I^I^^=observation point, 

N=number of time steps 
to be calculated) 


without affecting the numerical solution. When the initial field is im- 
pressed at z=0, the solution is straightforward and no special tech- 
nique is required at the sending end. 

The above discussion concentrates on the initial conditions 
prescribed in time at z=0 or normally incident from z 0. When the 
initial conditions are prescribed in space, i.e., 

E(z,0) = f(z) 

(62) 

g(z) 

then it is important to model the information at the time step At 
correctly. Carelessness here may sometimes lead to unexpected 
results. One example is shown in Fig. 5. Here f(z) is assumed to 
be a gaussian short pulse (see N=0) and g(z) is zero. It appears 
therefore that (62) for this case can be replaced by 
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However the numerical solution is quite strange as can be seen in 
the figure. One might attribute the numerical instability to the 
numerical solution itself, i,e., an unstable difference solution. 

Since we know that the numerical solution should produce an exact 
solution in this case as mentioned previously, there must be some 
other source of error. The truth is that the initial data e 1 in (63) 
is incorrect, because it implies the field at At is exactly^ identical 
to the fields at t =0. This is not physical. Consequently, a correct 
solution should not be expected. 

The correct modelling for this case is to employ the central 
finite difference[17, pp. 231] 



or 


E{z,-At) = E{z,At). 

Employing the difference equation (7) to At 

E(z,At) = - E(z,-At) + E(z+Az,0) + E(z-Az,0), 
v/e can eliminate E(z,-At) to obtain 

(65) E(z,At) = j [E(z+Az,0) + E(z-Az,0)] 

which is the corrected initial data at At for this case. In general, 
one should always ask, if the initial waveform is given at t=0, what 
would it be at At physically? By recourse to this kind of physical 
reasoning, unexpected solutions may be avoided. 

Finally it is always desirable and possible to incorporate the 
moving time window calculation technique to minimize unnecessary 
calculations in a difference solution. This can be carried out from 
a knowledge of the domain of dependence for the observation point. 

A computer program which is designed for Casey's problem[ll] is 
presented in Appendix A (Computer Program 1). Here the techniques 
discussed above are worked out in detail. 



CHAPTER III 

A DIFFERENCE EQUATION FOR STRATIFIED LOSSY DIELECTRICS 


In this chapter we shall extend the techniques developed 
earlier to the case of a stratified lossy dielectric medium. A 
difference equation appropriate for an interface grid point is 
derived, as done for the wave equation in the preceding chapter. 
The resulting difference equation is employed to solve three 
transient problems which were considered by other authors using 
standard transform methods. They include: a unit-step signal 

propagating in a homogeneous lossy dielectric, the transmission 
of a normally incident gaussian pulse from air to a conducting 
earth; and the reflected waveform of a sine-squared pulse from 
a three-layered lossy medium. In all cases normal incidence is 
assumed. 


It has been known for some time that the reflection of short 
pulses from stratified media offers a diagnostic tool in geophysics 
[25]. For example, Sivaprasad and Stotz[26] recently in a theoretical 
investigation concluded that the detection of water layers in dry 
earth at depths up to approximtely 10 meters is feasible. From 
another view point, the shielding property of the earth for electro- 
magnetic pulses has been of concern[27]. 

It will be shown that the solution of these transient problems 
can be easily obtained from a simple difference equation developed 
in this chapter. 

A. A Difference Equation 

It has been suggested in Section II.B that the lossy wave 
equation in a homogeneous medium[22] 


( 66 ) 





0 


can be approximated by the difference equation 


(67) 


ye 


e'?''’^-2E'?+e’?"^ 

_T 1 1 

At^ 


AZ^ 


ycr 


^n+l_En-l 
i ~ i 

2At 


= 0 


When the scheme 
(68) AZ = CAt 
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is selected, (67) takes the simple form 


(69) 




aAt 

2e 




which is then the difference equation for a homogeneous medium. 

To extend it to the case of a stratified medium, let us consider 
an interface boundary between two lossy dielectric media as shown in 
Fig. 6. Within each thin layer, the lossy wave equation is valid 
everywhere, hence it is true in the vicinity of the interface 
boundary. We write (66) at z~, z'*' respectively 






Fig. 6--Interface between two lossy dielectric media. 


At the interface the boundary conditions are 


E(z‘^,t) = E(z'.t), 



^2 9Z U-J 9z 


To derive a difference equation for the mathematical problem as 
posed in (70) and (71), we expand the fields at z'''+az 2 and z"-Az-| 
respectively by the Taylor series 
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( 72 ) 


E(z%AZ 2 ) = E(z‘".t) + 





E(z"-Az-|) 


E(z“,t) 



The second partial derivatives with respect to z in the above ex- 
pressions are first replaced by (70), the first partial derivatives 
with respect to z are eliminated by the boundary conditions in (71), 
and finally the second and first partial derivatives with respect 
to t are approximated by 
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AZ^ = ■ • At . 

\j ^ 2^2 

Note that li] , e] , a-i and u2» ^2* ^2 the constant values of the 

permeability, permittivity, conductivity of the thin layers to 

the left and to the right of the interface grid point i respectively.^ 

We see that (75) reduces to (69) when the electrical properties 
of the layers at both sides of an interface boundary are the same. 

It also simplifies to (52) when the inhomogeneous dielectric is loss- 
less, consequently, (75) can be used to obtain the transient response 
of a signal in an inhomogeneous lossy dielectric medium provided 
that the medium is approximated by many thin layers according to 

(76) Az. =-^At 
x/^i^- 

where the subscript i denotes the ith layer. In practice it is 
recommended that At be chosen such that aAt/2e <0.1 for accurate 
results. 

B. Numerical Examples 


Three numerical examples are presented to illustrate the 
applications of the difference equations in (75), / 

(1 ) Unit-step signal propagating in a homogeneous 
lossy medium 

We assume a unit-step signal is impressed at z=0, and it is 
desired to receive the signal at z=0.3, 0.6, and 1.05 meter in a 
medium with e = 16eo and o = 0.02 mhos/meter. With discrete time 
interval At = 1 nanosecond and Az = 0.075 meter, the waveform 
received at the corresponding observation point I = 5, 9, 15 are 
shown in Fig. 7. These results agree with those computed by Fuller 
and Wait[28] by analytical transform techniques. 

(2) Transmission of a qaussian pulse from the 
air into the earth 


This problem was originally considered by King and Harrison[27] 
in order to determine the shielding properties of the earth. They 
solved the problem by Fourier transform method. Here we obtain the 
time history of the pulse directly in the time domain. This approach, 
perhaps, is more efficient because a long excursion in the frequency 
domain has been avoided. Figure 8 shows the waveforms of the trans- 
mitted pulse at depths 0, 10, 25 meters below the surface of the earth 
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Fig. 8--Trans(ni tted part of a Gaussian pulse received at 
depth d in dry earth, tg is the signal front 
arrival time; t]=1.4l42 microsecond. 
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for the case of dry earth (e = 7€q. a = 0.001 mhos/m). The results 
obtained by King and Harrison are also shown. It is seen that the 
comparison of the results obtained by the two approaches is ex- 
cellent. 

(3) Reflection of a sine-squared pulse from a 
three-layered lossy dielectric medium 

We finally come to consider a three-layered stratified medium. 
Here a sine-squared pulse is normally incident on a layered earth 
(e = 7gq, a = 0.001 mhos/m) of thickness d with water (g = 81 gq, 
cf = 0.001) as the bed rock, and the reflected waveform of the 
diagnostic pulse is of interest. Figure 9 shows the results as 
obtained from a difference solution for several thicknesses of the 
earth layer. Our results again compare favorably with those obtained 
by Sivaprasad and Stotz[26] using the Fourier transform methods. It 
is interesting to mention in passing that the difference solution used 
about 1 second in computer CPU (central processing unit) time while 
the frequency domain approach requires about 20 seconds. This is 
not surprising because the latter requires responses at many fre- 
quencies to construct the time domain results, while the former in- 
volves only repetitive additions and multiplications, which are 
particularly suited for a digital computer. 

In summary, an explicit difference equation for application to a 
stratified lossy dielectric medium has been presented. Its validity 
has been demonstrated by three examples. A computer program which 
is designed for the last example is presented in Computer Program 2 
in Appendix A; little change is needed to produce the numerical 
r. suits for the other two examples. 


One final comment is that the electrical properties of dielectrics 
are always frequency-dependent. The propagation of a short pulse in a 
medium with constant values of permittivity, permeability and con- 
ductivity is thus a convenient rather than a realistic model. Never- 
theless, when the electrical properties of the medium do not vary 
significantly over the frequency spectrum of the pulse, then this 
model is a good approximation to the physical problem, which is often 
mathGfTiati Cd1 ly intractablG by transforrn niGthods* OnG can also omploy 
the difference equation to obtain the steady-state response for an 
inhomogeneous lossy dielectric medium by marching on in time long 
enough to allow the steady-state to be established. This has been 
illustrated in an example in the previous chapter. 
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CHAPTER IV 

A NUMERICAL SOLUTION TO SOMMERFELD'S "ABOUT THE 
PROPAGATION OF LIGHT IN DISPERSIVE MEDIA" 


A. Introduction 


One of the famous problems considered by Sommerfeld is a step- 
modulated sinusoidal signal propagating in a dispersive medium[l]. 

When Einstein in 1905 postulated that a signal can never be trans- 
mitted with a velocity greater than c, the speed of light in free 
space, a controversy was raised among the physicists at that time[2] 
because according to Stratton[22] , it was then generally believed 
that the group velocity was necessarily equivalent to the velocity of 
energy propagation; but it had also been known that in the absorption 
region the group velocity can be greater than c. This apparent con- 
tradiction was clarified by Sommerfeld who in 1914 used an elegant 
mathematical analysis to show that the front of the signal propagates 
exactly with c regardless of what group the velocity might be. Thus 
at a controversial time, this result conformed with one of the funda- 
mental postulates in Einstein's special theory of relativity. Sub- 
sequently, Brillouin[29] extended the analysis and gave a complete 
picture of the waveform of the precursor. He concluded that the main 
signal may be considered as arriving with the appropriate group 
velocity in the normally dispersive region, but that the group velocity 
loses its significance in the absorption region. 

In this chapter we shall present an extremely simple difference 
equation for thesolution of this historically famous problem. The 
difference equation is derived from an integro-differential equation 
for the electric field by a finite difference method. The solution is 
then generated in a recursive manner in time. Employing the difference 
equation, we have been able to simulate the propagation process in 
'ea]“time on a high-speed digital computer equipped with a CRT display 
facility. The dynamic behavior of the signal in a dispersive medium 
thus can be seen. In fact, we have observed vividly the formation 
of the precursor as the signal traverses the medium. We shall describe 
some of these observations and discuss briefly the numerical method. 
Extensive results will be presented to give a clear picture of the 
whole shape of the signal. 
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B. An Integro-differential Equation for the Electric Field 

Consider an x-polarized plane electromagnetic wave propagating 
in the z direction in a dispersive medium. The governing equations 
for the propagation process are the Maxwell's equations[30] 


(77) 


3E 

> 

at 


0 at 


3H 


az 


+ Ne|f=0, 


aE 


az 


^ = 0 


and the equation of motion for the electrons 


(78) 


ix 

dt2 


ds 

dt 


2, 

oj^S 

0 


^E 
m X 


where s represents the displacements of the electrons from their 
equilibrium positions. Here v, wq sre the electron collision 
frequency and the electron characteristic frequency respectively. 

Assuming the electrons are initially undisturbed, i.e., s, 
ds/dt = 0 at t = 0, the solution of (78) is given by [47] 


(79) s(t)=^^[ E^(z,B)sin h(t-S) exp[-p(t-3)]d6 

■*0 

/~2 2 

where h =Vuo"P p = Eliminating the magnetic field Hy in 

(77), the use of (79) yields an integro-differential equation for the 
electric field 


(80) 


a^E. 


at‘ 




az 


WpEx(z.t) = 


2“pP 


a(t-s)Ex(z,e)d& + 


where 



( h -P ) 


h 


‘'0 


b(t-e)Ex(z,B)de 


( 81 ) 


a(t-6) - cos h(t-e) exp[-p(t-$)], 
b(t-e) = sin h(t-g) exp[-p(t-6)] 
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^ 2 2 . 

and c = 1/v'yo^o’ “p “ />>ieo» c is the speed of light in free space, 

wp is the plasma frequency. The subscript x in E will be dropped 
hereafter. ^ 

C. The Difference Equation 


The numerical solution for this type of integro-differential 
equations has been discussed in Chapter II and is briefly treated 
here. The numerical method essentially replaces the second partial 
derivatives by the central finite differences 



E(z+A2,t)-2E(z,t)+E(2-Az.t) 


AZ 


2 


5 


(82) 


9 E _ E(z,t+At) - 2E(z,t) + E(z,t-At) 


ar 


At‘ 


To ensure a stable solution, it has been demonstrated in Section 
II.C that the electric field in the third term in (89) has to be 
approximated by the scheme 

(83) E(z,t) = [E{z,t) + |- E(z,t+At) + |- E(z,t-At)]. 


Finally, the integrals in the right hand side of (80) can be written 
as 


A(z,t) = 


B(z,t) 


= i- 

At 


-.t ^ 

•'t-At j 
/ft-« ft . 

\^0 Jt-At^ 


cos h(t-6) exp[-p(t-e)] E(z,e)dg, 


sin h(t-e) exp[-p(t-p)] E(z,e)dg . 


Evaluating the integration from t-At to t by the trapezoidal rule of 
integration yields the recursive relations 


(84) " g[KA(z, t-At )+SB(z, t-At )]+ E(z,t-At), 

B(z,t) =g[KB(z.t-At)-SA(z,t-At)] + ^ g KE(z,t-At) + ^ E(z,t) 
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where g=exp[-pAt], K=cos(hAt), S=sin(hAt). Combining the results of 
(82)-(84) and employing the notation = E(z,t) . , (80) is readily 

transformed into the difference equations t=nAt 


(85) 

where 


and 


cn+1 _ rn-1 . 1 

" "'^i ,.l 2 ‘-"i+l'"i-l 


1 - 


4 ^ 




aJ = 1 SE?‘\ 


B? = g[KB"''' - SAJ'^] + Y gKE""'' + 1 eJ 


T =tOpAt, = 2x^(pAt), T [(hAt)^-(pAt)^], 


Note that the scheme 


(86) z=cAt, c = 


has been used. In practice, it is recommended that the time increment 
At be chosen such that WpAt ± 0.1. 

When initial conditions are given, (85) can be used to obtain 
the time history of a signal of arbitrary waveform in the dispersive 
medium by marching on in time. 

D. Numerical Results 


1 


It is convenient to summarize Sommerfeld and Brillouin's results 
by considering a simulated experiment which we shall do shortly. A 
sinusoidal signal generator is located at Z=0; a detector of infinite 
sensitivity is set up a distance Z from the origin in the dispersive 
medium to detect the signal. For convenience we normalize the space 
and time variables by Z = wpz/c and T = wpt. At time T*0 the gen- 
erator is turned on. Then according to Sommerfeld and Brillouin[2] , 
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1, After a period of time T=Z (i.e., t=z/c), the front of the signal 
should be detected. The amplitude of the signal received im- 
mediately after is very small and can be described by 


(87) 


f(ZJ) = 


03 

2(T-Z) 

Wp, 

Z 


Jt(n/2Z{T-z)) u(T-Z) 


This is the first precursor and is valid only for small value 
of 2(T-Z)/Z. u(T-Z) is the unit-step function. 

2. At time 1=1 the second precursor arrives. 

r ^ 

3. At time Tg = Zc/vg {i.e., t = z/vg, Vg is the group velocity) 
the detector should detect an intensity equal to 1/4 the final 
intensity of the signal, which is true in the normally dispersive 
region. 

We shall attempt to verify these theoretical results by performing a 
numerical simulated experiment by the use of the difference equation. 

Instead of using a single detector, we set up four detectors at 
Z=0, 20, 30, and 40 respectively. The dispersive medium is assumed 
to be lossless and specifically = 1. Later, loss will be in- 
cluded. ° P 

Figure 10 shows the waveforms of the signal received at four 
detectors for the case u=0.2coo. It is seen that the signal is not 
distorted. This can be easily explained by considering the dis- 
persion relation (assuming eO“t) 

I 2 

(88) ^ ~ ^'1 ^ ^ ~ ^ “'^“o “ 

0) - (jJ ” ^ 

'"0 

which indicates that the medium is essentially non-dispersive for low 
frequency signals. The signal velocity thus can be defined clearly and 
is a little less than 0.7c. We can verify this velocity "experi- 
mentally" (in a simulation sense) by measuring the time T for the 
first peak of the signal and define without ambiguity the signal 
velocity 


V 


s 



c 


where Tq is the time when the generator reaches its first peak. Thus 
we have in this case 
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37^2 - 7.8 ^ 


0.683c, 

v^(Z=40) = 55 - 7,8 ^ 0.687c 


which agree with the theoretical result. In addition, one may also 
notice a minor disturbance in front of the main signal which is 
actually the precursor. It begins exactly at T=Z and thus confirms 
that the signal front propagates exactly with the speed of light in 
free space. We shall prove this later. 

Perhaps, a more clear picture of the precursor can be seen by 
increasing the frequency to ta=0.5a)o as shown in Fig, 11. Here a 
larger amplitude of the precursor is obtained as predicted by (87), 
but the signal no longer propagates without distortion. The leading 
portion of the signal loses its amplitude as it is traveling in 
the medium. However the waveform received at late times is un- 
distorted. A little arrow has also been added on the figure to 
indicate where the second precursor should occur as claimed by 
. Brillouin. However it cannot be identified as distinctively as the 
first precursor. 

The actual propagation process in space for this case is shown 
in fig. 12. As the signal front is traveling to the right in the 
medium, its leading edge breaks up Into little pulses. The deeper 
the signal travels, the more the numbers of pulses are formed, and 
the weaker is the amplitude on the leading peaks of the main signal. 
If a detector is located very far from the generator, one should 
receive more pulses in the precursor. Figure 13 shows this for a 
detector at Z=100. 

We now return to continue the experiment and increase the fre- 
quency to u) = oiQ. The received signals are shown in Fig. 14. Here 
the waveform is distorted to such a degree that it has nothing in 
common with the original signal as it was launched (i.e., at Z=0). 
Consequently, there will not be a suitable definition for the 
velocity of propagation in this case. But one thing is clear: the 

front of the signal always arrives with a speed equal to c. 
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Fig. 11 --Received waveform at Z of a step-modulated signal in a 
dispersive medium. of.Buiq, »o=“p> T=0.1N. Arrow 
indicates normal position of the second precursor. 
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Fig. 13— Received waveform at Z=100 of a step-modulated 
sinusoidal signal in a dispersive medium. 
f T — O.IN. 
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We further increase the frequency to o)=1.5u)o. The waveforms 
received at each detector are shown in Fig. 15. In this case the 
signal is building up gradually and there is no clear distinction 
between the precursor and the main signal. One useful characteri- 
zation, the build-up time tg = z/Vg, can be defined however. 

According to Brillouin, it should occur at a time when the de- 
tector detects an intensity equal to 1/4 the final intensity. Notice 
that this also applies to previous cases, when the signal is built up 
gradually, such as the case u=0.5a)^ in Fig. 13. 

Finally, we consider a specific loss p/a>o=0.1 in the medium and 
repeat the experiment. The received signals are shown in Figs. 16- 
18 for (d/uq = 0.5, 1.0 respectively. These results indicate 

that the presence of loss simply attenuates the amplitude of the 
signal. This is particularly apparent in the case w=0.5a) by com- 
paring Fig. 17 and Fig. 11. 

We present here a simple proof, using the difference equation, of 
the fact that the front of a signal always propagates with velocity 
c and that the front also preserves its amplitude in the course of 
propagation. Let us consider the space-time grid pattern as shown 
in Fig. 19. Here the rectangular net is constructed such that 
Az=cAt (z=iAz, t=nAt). The signal generator is located in the first 
column i=0. Now we apply the difference equation (85) to all grid 
points, along the diagonal line AB and determine 


(89) 



1 + 


(t»pAt)' 


E° 

0 


(n-1 ,2 ,3, • • • ; i-n) 


which describes the course of propagation for the signal front E 
Equation (89) can be approximated by 


0 

0 * 


e; - [1 - n(«pSt)2] e". 

= 0 - T ‘ E? 

where t is the time of flight for the signal front from the- source 
(Z=0) to the observation point and hence is finite (actually t = 
z/c as we shall see). We then have in the limit At ->■ 0 that 

(90) = Ep (n=l,2,3,...; i=n). 
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Suppose n=l , we have e]=Eq. The interpretation is that the 
signal front E®, after spending At seconds, has travelled a distance 
Az meters to e] Without losing its amplitude. Hence the velocity of 
propagation for the front is 


( 91 ) -s.f. 

Thus the signal front not only propagates with c, but also preserves 
its amplitude in the course of propagation. For example, if an 
ideal unit-step signal of zero risetime is suddenly applied at z=0, 
then at any point in a dispersive medium one should detect an electric 
field of 1 volt/m arriving exactly with the speed of light in free 
space. In this connection it is found that the statements in some 
books such as "••• so that the process starts always from zero ampli- 
tude."[22, pp. 337] and "••• However, the amplitude of the first 
impulse is zero * •♦"[31] is true only when E°=0. This is certainly 
the case for a step-modulated sinusoidal signal considered by 
Sommerfeld. 

E. Summary 

In summary, we have found "experimentally" that the front of the 
signal indeed always propagates with the speed of light in free space. 

In addition, the front of a signal preserves its amplitude in the 
cor'^se of propagation. The physical reason for this is that, for the 
Lorantz model dielectric, the relative complex dielectric constant 
tends to unity for high frequencies. As the signal travels through 
the medium, its leading edge breaks up into little pulses. These pulses 
are picked up first by the detector and are known as the precursor. 

When 0 ) « (Dp, its amplitude is small and the signal propagates es- 
sentially without distortion. As (d increases, so does the amplitude 
of the precursor. But now the leading portion of the main signal 
suffers a decrease in amplitude as it traverses the medium. When 
(D (Dq, the signal is severely distorted. Hence the concept of 
velocity of propagation is not useful here even though the front always 
travels at speed c. When (d » (Dq, the signal received at a detector 
is built up gradually and there is no clear distinction between the 
precursor and the main signal. A build-up time can be defined however 
and it is related to the group velocity. The presence of loss in the 
medium simply attenuates the strength of the signal. 

A computer program is presented (Computer Program 3) in Appendix 
A which can be used to generate the results in this chapter. 
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CHAPTER IV 

TRANSIENT SIGNAL PROPAGATION IN COLD STRATIFIED PLASMAS 


A. Introduction 

The propagation of a transient signal in cold plasmas has been 
investigated extensively in recent years[e.g., 4, 32, 33, 34], To 
analyze this class of problems, Fourier and Laplace transform methods 
are usually employed. Although these standard methods apply to a 
number of cases, their failure to accormnodate others has warranted 
additional and new methods of solution. For instance, Bowhill[6] 
suggested the multiple-scattering technique while Field[7] proposed 
the method of characteristics. 

This chapter extends the previous one to a stratified plasma 
where a constant magnetic field parallel to the direction of propa- 
gation is also assumed. By approximating the plasma medium with a 
large number of thin layers, it is shown that the time history of a 
signal of arbitrary waveform in an inhomogeneous, anisotropic {because 
of the constant magnetic field), lossy, cold plasma can be easily 
obtained from difference equations by marching on in time. Only the 
case of normal incidence is considered. 


In Section B the coupled integro-differential equations for the 
electric fields are first derived from the Maxwell's equations 
together with the equation of motion for the electrons. These 
equations are next transformed into difference equations in Section C. 
Numerical techniques and examples are presented to demonstrate the 
versatile applicability of the resulting difference equations. 


Coupled Integro-differential Equations for the Electric Fields 


Consider a plane electromagnetic wave propagation in the z 
direction in a homogeneous plasma in which a constant magnetic field 
Bq parallel to the direction of propagation, is assumed. The 
governing equations for this problem have been derived from Maxwell's 
equations by Stratton[22, pp. 328] and can be written in complex 
form 


(92) 

with the corresponding equation of motion for the electrons 
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(93) (v + 


where E = Ex+iEw, H = Hx+iHy, v = dx/dt + i dy/dt, and wq = eBo/m. 
Here N, e, m, v", uc are the electron density, the charge of an 
electron, the mass of an electron, the electron collision fre- 
quency, and the electron cyclotron frequency respectively. 

Assuming the electrons are initially rest, the solution of 
(93) is given by 


(94) 


v(z,t) = f 


E(z,3)exp[-(v+iu )(t-6)]de. 


Eliminating the magnetic field H in (92), the use of (94) yields 
the integro-differential equation 


(95) — 

at 


■^+ u^E(z,t) = a)^(v+ia> ) E(z,6)exp[-(v+iio )(t-e)]de 
P P c Jo c 


2 2 

where c = 1 /J po^o » <^p " /"leo and c is the speed of light in free 
space, top is the plasma frequency. Equation (95) is further de- 
composed into real and imaginary parts to obtain the coupled time- 
dependent integro-differential equations for the electric fields 


3^E p 9^E p 2 

(96a) P - — ^ + (o^E (z,t) = (0 

at^ az"^ P 


P J 


a(t-e)E^(z,6)dg 


- o)p I b(t-B)Ey(z,e)d6, 


(96b) 


a^E.. o a^E 


ar 


2l^..!Ejz,t) = 


rt 


az 


2 ’ fjj 

p y 


(jS 


P J 


a(t-B)Ey(z,B)de 


+ to' 


p j 


b(t-B)E^(z,6)dB 


where 



a(t-fj)-[vcosu) (t-0)+w si nto (t“B)]exp[-v(t-$)] , 
b(t-6)=[u)^cosa)^(t-6)-vsinu)^(t-B)]exp[-v{t“p)]. 
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Here we can see clearly the coupling of the electric field com- 
ponents which occur in the last terms in (96). 

To visualize this coupling, let us consider an x-polarized 
plane electromagnetic wave in free space normally incident on a 
plasma medium. Before the signal front reaches the air/plasma 
interface, the electric fields are, of course, governed by the 
wave equation, i.e., by setting wp=0 in (96). Since there is no 
Ey component in the incident wave, Ey is identically zero in this 
period of time. At the instant Ex enters the plasma, an Ey field 
is set up by the last term in (96b). This induced electric field 
then propagates both into the plasma and back to the free space. 
Therefore, the reflected wave has an Ey component. It is clear 
that the last terms in (96) are responsible for converting the 
polarization of the incident signal. When the kernal b(t-g) is 
zero, i.e., there exists no constant magnetic field, then (96) 
becomes 

2 2 t 

(97) ~ ” V [ exp[-v(t-e)E(z,3)dS 

sr 9z^ P P Jo 

where E can be x or y polarized and the subscript is dropped. Here, 
the polarization would be preserved since the equations are now 
u. coupled. One can also obtain this integro-differential equation 
fi-om (80) by setting the characteristic frequency wo=0. 

C. Difference Equations for the Electric Fields 

To transform (96) and (97) into difference equations, the plasma 
medium is first divided into a number of thin layers of equal thick- 
ness 

(98) Az=cAt 

in which the time increment At is arbitrary but small; its choice 
will be clear later. We then attach to each interface boundary a 
grid point as shown in Fig. 20. 
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(1 ) Homogeneous plasma 

A procedure to derive the difference equations for this case 
has been described in Chapter II. The numerical method essentially 
replaces the second partial derivatives by the central finite dif- 
ferences 


(99) 


^ _ E^(z.t+At)-2E^(z,t)+E^(z,t-At) 


9t‘ 


At" 


where c stands for both x and y. To ensure a stable numerical 
solution the electric fields in the third term in (96) are ap- 
proximated by three consecutive times 

(100) E^(z,t) = 1 [E^(z,t) + i E^(z,t+At) + 1 E^(z,t-At)]. 

Finally, the integrals in the right-hand sides of (96), (97) are 
wri tten as x / > \ / 


{ rt~At rt \ 

Ar(z,t) = + 

'■'o •'t-At / 


a(t-6)E^(z,e)d6, 


It At) 


D(z, 


" At ( f +f ] exp[-v(t-8)]E(z,B)d$. 

v-^0 H“At/ 


Evaluating the integration from t-At to t by the trapezoidal rule 
of integration yields the recursive relations 
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A"i = g(KA"^^ + j [aE". + (0S+aK)gE":^]. 

(101) b". = g(KB"'^ - SA"^^) + i [0E". + (eK-aS)gE"‘h 

nH - X 1 rc1 X ^r^-li 


where K=cos0, S=sin6, g=exp{-a), a=vAt, 0=uicAt and the notation 
E^^ = has been used. 

Z=i AZ 


Combining the results of (98)-(101), Eqs. (96) are readily 
transformed into the explicit difference equations 


En+1 ^ _ ^n-i ^ t^[e” + e>i . iZ^n. + t^(a".-b".)], 

XI XI '■xi +1 X 1-1 2 XI ' XI yr-*’ 


-n-1 


yv 


(102) 


2 

rn'f'l _ cO“l 1 riTF^ X IE X ^ xd^ M 

i+1 ^ ^y i-l “ T ^yi ^ " (\i^^xi^^ 


where t = WpAt 


3 


n 


1 



and similarly (97) is transformed to the difference equation 


(103) 


:n+l _ _^ri-l 


n-1 X v^r^n x t x 2 -.n^ 


(2) Stratified Medium 

When a grid point is located at the interface between two dif- 
ferent plasmas, or the air/plasma interface, strictly speaking 
Eq. (102) no longer applies. Not only is the finite difference in 
(99) for the z variable no longer valid, but even the integro- 
differential equations themselves are undefined here. Fortunately, 
appropriate difference equations can be derived. 

Without loss of generality let us consider the interface 
between layer 1 and 2 in Fig. 20. Within each thin layer the 
integro-differential equations are continuous everywhere, hence 



they are valid in the vicinity of the interface boundaary. Thus 
we write (96) at z"*", z~ respectively 


n “ C ' “t" QJ 


(104) 


sr 


9 Z 




9^E (z'.t) 9^E.(z',t) . 

-C + oj 


9t‘ 


-^2 *“p, E 5 U .t) = »;^F,[v,.E^(z-,t)] 


where the subscripts 1 and 2 denote the thin layers to the left and 
to the right of an interface boundary respectively. Each of the 
above expressions actually stands for two equations; the functionals 
F] , F 2 are understood as representing the source term in the right- 
hand sides of (96). 






Fig. 20" Interface between two 
lossy cold plasmas. 

At the interface the boundary conditions are 
E^(z'^,t) = E^(z".t), 


( 105 ) 
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To derive the difference equations for the mathematical problem 
as posed in (104) and (105), we expand the fields at the points z"^+Az 
and z“-Az by the Taylor series expansions 


+ . 8E (z''’,t) 

E^(z +Az,t) = E^(z ,t) + ^ — Az + 




3Z 


2 , 


(106) 


9E (z‘,t) 

E^(z -Az,t) =E^(z ,t) ^ AZ + 




3Z 


in which the higher order terms have been neglected. The second 
partial derivatives in these two expressions are first substituted 
by (104). The first derivatives are next eliminated by the use of 
the boundary conditions in (105). Finally we replace the second 
partial derivatives with respect to t in the resulting equation by 
the finite difference in (99). Performing all these steps, we 
obtain the difference equation 


(107) E^(z^,t+At) = -E^(z^,t-At)+E^(z^+Az,t)+E^(z'-AZ,t) 

2 + T^) E|(z^,t) 


+ ^{T^F^[v^.E^(z",t)]+T^F 2 [v 2 .E (z'',t)]} 


where T]=ci3piAt, T 2 =wp 2 At. A is added to indicate that this term 
must be approximated by the scheme in (100) for numerical stability 
considerations. 


The above results indicate that for an interface grid point one 
simply evaluates the lower order terms (i.e., the plasma and the 
source terms) on both sides of the interface and then take the 
average value of them. Therefore, (107) can be written explicitly 
as 


:n+l 

'xi 


-E 


n-1 

xi 


n[E'^..,+ e" 

Xl+1 X 


i-1 


-yE 


XI 


^ j=l 


j' xji yjv-“ 


e"!’ = -e":'. 
yi yi 




+e*^ 
i+1 y 


yE". + 


1 


i-1 ’"yi 


2 

2 I 

^ :=i 


y yji xj!^-! 


(108) 
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where 


and 


''"jl " ^ (es+c.K)g E^:'], 

\J 

T, = u At, T„= (D At, a. = v.At, g^=exp(-a.), 0 = w^At, 

* p-| ^ ^2 J J J J ^ 

K = cos e,s = sin e, y= j (t^+t,), n = — — ^ — 

^ + jy ^ 

Similarly, we can write down the difference equations for an isotropic 
lossy plasma from (97) 


(109) 


Eii+1 , . Erl 

1 1 


* nCEj,! + E;_, 


YE*^ + — 

''i 2 



2 



It is seen that (108) reduces to (102) when the plasmas at both 
sides of a grid point is the same. Furthermore, since (108) is derived 
from a canonical problem of an interface boundary, it applies to other 
interface and hence a multitude of boundaries as well. Consequently, 
(103) can be used to obtain the time history of a signal in an in- 
homogeneous plasma provided the medium is approximated in a step- 
wise fashion. 

D. Some Practical Aspects of the Difference Solution 

To illustrate some pertinent techniques for the application of 
(108), let us consider the space-time grid pattern as shown in 
Fig. 21. Here the space is divided into thin layers of thickness az. 
The time is quantized to an increment At (recall that Az*cAt). The ’ 
field at a grid point is denoted by EO. Distance and time are 
defined as z=(i-l)AZ, t=nAt. 

For the time being we assume the incident signal is i mpressed 
at z=0 (i=l). It is desired to obtain the transient response up to' 
a time T=a)nt at the observation point Z=[opz/c. For simplicity we 
consider that the plasma is homogeneous, lossless and isotropic, 
although the discussions that follow apply to a general situation 
as wel 1 . 



62 


i = 1 2 3 4 ... I 



C 


Fig. 21— Space-time grid point pattern, 
z=(i-l)Az, t=nAt. 


The difference equation for this case is the simplest one and 
is given by 


+- 

1 


^•" i+1 1-1 


1 + 4 T 


{i=2,3,4,-..,I-l; n=2,3,4,-.-,N) 


where E may be x or y polarized. At the column i=l we have the initial 
condition 


(111) E" = f(nAt) (n=0,l,2,-..,N) 

where f(t) is the source signal impressed at z=0. Since we have 
assumed the plasma is initially undisturbed, we have at n=0,l 
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E? = 0 (i > 2). 
e[= 0 (i>3). 

I 

Now the field at E2 is no longer zero because the signal front at 
E^ has reached this point. We can however apply (110) and determine 

(112) eJ = ^^1- 



Here the fact E~p=0 has been used. This follows from causality of 
the fields. ^ 

We finally select an appropriate value of r, and employ (110) 
to generate the solution up to the time step 

(113) 


At each step the response is simply picked up at the observation point 


(IM) 


^obs Az T 


+ 1 . 


The far-end boundary can be ended at 
I = iobs + ^ (N-iobs> * ' 


without affecting the numerical solution. This follows from the 
domain of dependence for the observation point. For this very same 
reason, only those grid points inside and on the trapezoid ABCD in 
Fig. 21 need to be calculated. 
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The condition in (114) is good for an infinite plasma medium. 
When the plasma is bounded by a perfectly-conducting boundary, one 
can simply set 


(116) eJ = 0 (n=0,l,2,-..N) 


at the end point of the conducting boundary. 

When the plasma slab is bounded by the free space, then the 
assignment 


(117) e{ = ej:] 

would simulate the infinite free space. 

There remains an important case for this class of problems; 
that for which the incident signal , rather than being impressed at 
z=0, is incident from z<0 in the air. It is then numerically ef- 
ficient to assume the signal front arrives at z=0 (i=l) at the 
instant t=0. This is equivalent to assigning again the initial 
fields to column i=l , although a different physical problem is 
encountered here. 

Now the air/plasma interface can be conveniently chosen at 
i=2. While (110) remains valid for all interior points in the 
plasma, the difference equation for i=2 is 


(118) . -E 5-2 . (E"-’ * S - 1 E^-') 

’ ^ 8 ^ 

where S is literally equal to E. \ but has a different physical 
signifii:ance in this case. Physically the quantity S means the 
total field at a point where it is in the air. Thus it has two 
components: the incident field and the reflected field. Refer to 

Fig. 21^ and suppose (say n=3) is being evaluated; then the 
total field in the air is given by 


(119) 


S = Primary field H + Reflected field from G 





). 


Here eP ^ is the incident field as discussed earlier, E^"^ is the 
-total field at G and E*j'"-^ is the incident field to 6 from A. 
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Hence the reflected field from G is equal to the difference be- 
tween these two quantities. The reflected signal is then simply 
picked up at i=2 according to 


(120) E'^(nAt) = e" - (1 < n < N) 

The problem is then completely solved. 

We now want to comment on the selection of the time step size 
At. In practice, it is recommended that At be chosen such that 


(121) (u At, w At, vAt) <_0.1 
P ^ 

whenever possible. The reason for the choice (121) is solely due 
to numerical stability and accuracy considerations. One should recall 
that we have obtained the recursive relations for the integrals by 
trapezoidal rule of integration as in (101). When the step size is 
too large, it is not clear whether the error introduced at one step 
will jeopardize the stability of the solution at later times. Linder 
such circumstances, we have had recourse to intuition and numerical 
evidence, since we have not derived a stability condition for this 
type of equation. 

E. Numerical Examples 

The validity of the resulting difference equation is demonstrated 
by several examples selected from the literature. 

The first example is based on the work of McIntosh and 
El-Khany[35] in which a group velocity concept was used to synthesize 
a "chirp pulse", and the transform methods were applied to obtain the 
actual "compressed" pulse. Starting with their signal waveform, 
we have calculated the waveform received at Z=15 in a homogeneous 
isotropic plasma for two cases v/wp = 0, 0.1 (see Fig. 22). Although 
only the lossless case was considered in [35], we see that pulse 
enhancement still occurs in a lossy plasma. More recently, pulse 
compression in bounded dispersive media has also been demonstrated 
by transform methods[36]. With the aid of the difference equation, 
the response of the signal constructed from a number of simple models 
can be easily computed in various plasma conditions. This practice may 
lead to computer-aided design procedure for pulse compression studies. 




Fig. 22— Pulse enhancement at a receiver 
Z=15, T=NwpAt, o)pAt=0.2 
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The next example compares the difference equation solutions for 
the reflection of a unit-step signal with results obtained by Wait[37] 
(for the lossless case only) by transform methods. Figure 23 shows 
the reflection of a unit-step signal normally incident on a homogeneous 
plasma half-space. The curves are presented for v/iop = 0.0, 0.1, 0.5 
and 1.0. In an accompaning table the difference solution for two 
different time step sizes ( = 0.1, 0.01) is compared with the 

analytic solution[37] for the lossless plasma case. The convergence 
of the numerical solution is clearly indicated. As was mentioned 
earlier, if the electric field in the plasma term is not weighted 
at three consecutive time steps as in (100), then the difference 
solution is unstable. This solution is given in parenthesis, and is 
unrecognizable after T > 5. 

The third example is to calculate the waveforms of a short 
gaussian pulse reflected from an inhomogeneous plasma having a 
linear electron density profile ((o 5 = 2 xl 0 l 0 z) and constant collision 
frequency with a vertical magnetic^field. This is shown in Fig. 24. 
This problem was originally considered by Hill and Wait[38] using 
transform methods, and numerical procedures were used to obtain the 
results. The difference solution are obtained by the use of (108) 
by marching on in time. Comparing Fig. 24 with the corresponding 
results in [38], we conclude that the accuracy of the difference 
solution is indeed excellent. 

The fourth example (Fig. 25) shows the reflection of a unit-step 
signal normally incident on a lossless plasma slab for 4 different 
normalized slab thickness Z = 30, 5, 2, 1. Notice that the first 
case is essentially the half-space case. Wait[37] has obtained an 
analytic expression for this problem in the form of an infinite 
srries involving Bessel functions. In the presentations of his 
curves, there are still noticeable oscillations in the reflected 
waveform at T > 10 for the cases Z = 1, 2. However, our difference 
solution does not exhibit these oscillations. Antonucci[45] had 
constructed an artificial TEM line to simulate the plasma and also 
noted that his result did not agree with Wait in these cases because 
his results showed no oscillations at late times. 

The fifth example shows the response at Z = 100 of a step 
carrier signal at Z = 0 in a homogeneous lossless cold plasma. Notice 
that the signal builds up gradually and oscillates about the steady 
state values at late times. This phenomenon is predicted analyti- 
cally[32] and confirmed experimental ly[48, 49] in related problems. 

For our final example we show the spatial variation of a step 
carrier at z=0 for three frequencies w/un = 1, 1.2, 1.4 propagating 
in a lossless plasma at several normalized times (see Fig. 27). When 
0 ) = 0 ) (and to<o) ), we see that the waveform is strongly attenuated as 



it propagates to the right. But the wavefront of the signal always 
propagates with the speed of light. Other velocities become dif- 
ficult to define. When u > wp, the signal builds up gradually 
depending on the frequency of the propagating signal 

The above results can be obtained from two computer programs 
presented in Appendix A. Computer Program 4 calculates the propa- 
gation and reflection of a unit-step, step-modulated sinusoidal 
signal, and phase-varying rectangular pulse[35] for a homogeneous, 
lossy, half-space or slab geometry. Computer Program 5 calculates 
the reflected waveforms of a unit-step, step-modulated sinusoidal 
wave, and a short gaussian pulse from an inhomogeneous plasma[38] 
or homogeneous plasma half-space with a vertical magnetic field. 
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Fig. 23--Reflection of a unit-step signal from an isotropic 

plasma half -space versus N. T=No)pAt, d)pAt=0.1. 

{Entries in parenthesis show the case in which 

E is not approximated by 3 successive values) 

X 
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Fig. 24- -Reflected waveforms of a gaussian short pulse normally 
incident on an inhomogeneous plasma having a linear electron 
density profile with a vertical magnetic field versus N, 
t=NAt, At=17.4 nanoseconds. 

(a) v=10^ sec"l , and w^=7xl0® for 

(b) v=10^ sec"l , 

(c) v= 2 x 106 , 

(d) v=0.5x106.e1"c=o,eJ^c= — - - - exp[-(N-ll )^/2k2] 

'^2tt(0.02I<)2 
{k=4. 0<N^22). 
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Fig. 25— Reflected waveform of a unit-step signal normally 
incident on a homogeneous cold plasma slab of thickness Z. 

T=0.1N ; 
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0.0 l.O -l.O 0.0 1. 



Fig. 26— Received waveform at Z=100 for a step-modulated sinusoidal signal 
in a lossless isotropic plasma, T=0.1N. 


ro 


Fig. 27--Spatial variation of a sinusoidal signal in a 
lossless plasma at several normalized times. 

T=N(d At, 0) At=0.2 



CHAPTER VI 

UNIT-STEP SIGNAL REFLECTED FROM A DEBYE DIELECTRIC 


A. Introduction 


Fellner-Feldegg[40,41] and others[42] have recently demonstrated 
quantitatively the feasibility of determining the dielectric constant 
of dielectrics over a wide band of frequencies by utilizing time- 
domain ref 1 ectometry , i.e., by measuring the reflection of a unit- 
step signal from a sample of the dielectric terminating a coaxial 
transmission line. In order to investigate the possibility of a 
direct determination of the low frequency permittivity, the high 
frequency permittivity, and the relaxation time spectrum of a 
dielectric in the time domain, theoretical calculations of the 
reflected waveform of a signal must be made. This has been done 
in the past using Laplace transform method[43,44]. In this chapter 
it is shown that this problem can also be solved effectively by the 
finite difference method. An explicit difference equation is derived 
from the governing time-dependent integro-differential equation for 
the electric field in a Debye dielectric with a single relaxation 
time. Employing this difference equation the time history of a signal 
of arbitrary waveform can be generated by marching on in time. Thus 
the response of a nonideal unit-step signal with finite rise-time 
can be computed. This may give a better comparison wih the meas- 
urements, because in practice a signal always has a finite rise-time. 
As an example of the treatment of the Debye dielectric case, we shall 
also calculate the reflected waveform of a unit-step signal incident 
on an ice layer on water. Here the water is assumed to be a Debye 
dielectric. 

B. An Inteqro-Differential Equation for the Electric Field 
Consider the Maxwell's equations 




9H 


at 


5 

az 


= 0 


which is equivalent to 
(123) 


2 

a'^D. 


a^E 

1 ^ “^x 


at 


Vr 


az 
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For a Debye dielectric with a single relaxation time, we have[12] 

(124) D^(z,o>) = e*(o.)E^(z,(o) 

where 

(125) e*(u)) = e + ^ ^ 

1 + ju)T^ 

and E£,eoo are the low and high frequency permittivities, is the 
relaxation time. Then (124) and (125) are transformed into the time 
domain 


Dj^(z»t) = 


.OD 

J — oa 


e(t-B)E^(z,B)dB 


9 


e(t) = 


and thus. 


e <s(t) + 




(126) D^(z.t) = + 


£ GO 


exp 


f 1 

- E fz,e)u(t- 


(t-B)dg. 


Here the lower integration limit can be changed to 0 because of 
cau'^ality of the fields, and the upper integration limit can be 
c'. anged to t by omitting the unit-step function in the integrand. 

Ije differentiate (126) twice with respect to t and obtain from (123) 


(127) 




E^(z,t) 


where 



For completeness the static conductivity oq of the Debye dielectric 
(as in the case of sea water) is included. (127) governs the 
electric field component of a plane electromagnetic wave in a Debye 
dispersive medium. 
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C. The Difference Equation 

In this section we shall derive the difference equation from 
which one can obtain the reflected waveforms of a signal when it 
is suddenly applied to the air dielectric interface. The geometry 
of the problem is shown in Fig. 28. The Debye dielectric is first 
approximated by thin layers of thickness 


( 128 ) At 

/ y^e 
*>/ 0 “ 

where At is arbitrary but small. There are two difference equations 
needed: one for the interior point in the dielectric, and one for 
the interface point. However no calculation is needed for the points 
in the air, as discussed in Section II. H, for the signal incident 
from z < 0. 


Employing the notation eJ = E(z,t)^^^.^^ 

t=nAt 


the difference equation 

for interior points can be written directly from (127) 

^n~l 

( 129 ) ~ J i - c 


2 - 2e".e;_^ 


At 


Az; 


-n+l ^n-l 




Here sj denotes S(z,t)^^.^^ which is defined as 


t=nAt 


/ V 1 . 1 /ft-At 

S(z,t) = exp[-o)^(t-B)]E(z,B)dB = ^ 


f ) 

•*t“At/ 


Thus , 
(130) 


• exp[-Wjj(t-B)]E(z,e)dB 


sj 


gsj ^ +1 (eJ + gE?‘^), g = exp(-w^At). 


Employing the scheme (128), Eq. (129) can be simplified to 
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+ a(w^At)^E9 -a(w^At)^s'j' 

where o=a +ao)^E„. Equation (131) applies to all grid points in the 
Debye dielectric. 

The difference equation for the interface point can be derived 
by the same procedure as used several times previously. Refer to 
Fig. 28. In the vicinity of the air/dielectric interface we have 


(132) 


S^E(z~.t) __2 9^E(z',t) _ 

2 ^ 2 ~ 
3t^ 



_J 



. f- -aA(z^t).F(z^t)=0 

8t 3Z » ^ ° 


where 

F(z'^,t) = aw^ 


>t 

•'o 


-%(t-e) + 

e E(z ,g)dp 


AIR 

Z“ 

DEBYE 

Z+ 

: DIELECTRIC 

% 

w ^ 

-AZ,- 




Fig. 28 — Interface between air-Debye dielectric. 


The boundary conditions at the interface are 



E(z‘",t) = E(z‘,t), 
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(133) 

3 E(z^«t) _ 3E(z~.t) 
3Z 3z 


To derive the difference equation from (132) and (133), we ap- 
proximate the air and dielectric by thin layers with different thick- 
ness in each region such that 

Az-, = C At, 

I 

(134) 


AZo = C At. 
2 “ 


We then expand the field at z +AZ 2 and z~-Az^ by the Taylor's series 
E(.%.z„t) = E(z\t) . .z, * i . 


035) 


3Z 


E(z--.z,,t) = E(z-.t) - AZ, _ 


The second partial derivatives in the above expressions are first 
substituted by (132), the first partial derivatives are eliminated 
bv the use of the boundary conditions in (133) and second partial 
derivatives with respect to t are approximated by the central finite 
differences 


3^E(z-,t) ^ E(z~.t+At)-2E(z-.t) + E(z~,t-At) 


3t 


At" 


Completing all these steps we obtain the difference equation for the 
interface point 


(136) aEj'^^ 






where 
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Since p36) will be used only on the air/dielectric interface at the 
point i=2, we can write it specifically 


(137) 






+ f [a(wQAt)^E2-a(to^At)^S2]. 


Here the quantity e" represents the total field for a point in the 
air. Thus, for the unit-step signal 

Total field e!^ = Incident field + Reflected field 

= 1.0 + (e""^-1.0) 


The complete algorithm for the solution of a unit-step signal 
reflected from a Debye dielectric is then given by 


n=0 

n=l 


( 138 ) 


n=2,3,4--' 


rn _ 1.0 (i=l) 

“"i ■ 0.0 (i>l) 


i”. - 


1.0 ( 1 = 1 ) 

CO 

0.0 (i>2). 




037 ) (i=2) 
(131) (i>2) 
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D. Numerical Results 


Figure 29 shows the reflected waveform of a unit-step signal 
from Debye dielectric with e£=80, eoo= 8, and aoTo/eo=0» 11.3, 113, 

1130. These results agree excellently with those obtained by employing 
the Laplace transform method[44]. Note that the amplitude of the 
waveform for early and late times corresponds to the reflection coeffi- 
cient for high and low frequencies respectively for the lossless 
case[44]. 

We next consider an ice layer on water as shown in Fig. 30. The 
thickness of the ice is 30 cm and its permittivity is ei=3. The water 
is assumed to be a Debye dielectric with ej^=4.5, Eoo= 81 and is the 
relaxation time. A unit-step signal is suddenly applied to the air/ice 
interface; after a time of 

j _ 0^30 _ X lo"^ seconds 

J3 

the front of the transmitted signal reaches the ice/water interface; 
transmission and reflection occurs, the reflected signal takes another 
Tq seconds to return to the air/ice boundary. Thus in the period 
0<t<2To the waveform received in the air is solely due to the re- 
flection from the air/ice interface. Now the reflected signal from 
the ice/water will be reflected again by the air/ice boundary, but 
this signal will not arrive the ice/water interface until t=3To. 

Thu*" we may neglect the multiple reflections in the ice by terminating 
tl.3 transient response at 3 Tq. The result of such a calculation is 
shown in Fig. 30. This problem could not be solved practically 
without the techniques discussed in Section II.H because then a large 
number of grid points (about 30000) is needed to model the ice. This 
is because of the exceedingly small relaxation time of the water 
relative to the time for propagation through the ice. 

A computer program for this problem is given in Computer Program 
6 in Appendix A. In this program the static conductivity oq is also 
included. The result in an earlier example (Fig. 29) can also be 
obtained from this program by setting the permittivity of the second 
layer, i.e., the ice, equal to 1. 
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REFLECTED STEP RESPONSE 



Fig, 30 — Reflected waveform of a unit-step signal normally incident 
on an ice layer of thickness 30 cm on water. tQ=3.46 nanoseconds. 




CHAPTER VII 

DISCUSSION AND CONCLUSIONS 


In the previous chapters we have shown how finite difference 
methods can be applied to linear dispersive waves. But they are 
even more powerful for nonlinear waves, which unlike the linear 
case, cannot be handled by transform methods. One of the nonlinear 
dispersive waves which caught the attention of many investigators 
in recent years is the Korteweg-deVri es equation[13, 46] 

(139) 

3X 

This equation was originally suggested for shallow water waves (1895), 
but it has been found to describe the long-time evolution of small, 
but finite amplitude nonlinear dispersive waves, such as the magneto- 
hydrodynamic waves and ion-acoustic waves in cold plasmas. 

In 1965, Zabusky and Kruskal[13] proposed the difference equation 
for (139) 


uf’-u?-’ 

2At 


+ ± n^ n x ^i-H~^i-l . .2 1 

+ 3 (u.^.^+u.+u._^) ^ + 6 ^ 


1 


(u!. 


i+2 ^'^i+l'''^^i-r’^i-2^ "" ^ 


wrn'ch can be written as 


(140) u!.' • = uv 


n+1 _ ,,n-l 1 At 




2 

6 At 


f - 2-HI * - < 2 >- 


Ax 

Given the initial conditions, (140) can be used to calculate the 
evolution of the initial waveform. Such a calculation for the case 
of the periodic initial condition 


(141 ) u(x, t=0) = COSirX 
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is shown in Fig. 31, given by Zabusky and Kruskal . We see that the 
initial waveform breaks up into eight pulses eventually. After that 
time thpe pulses, now called solitons, begin to move without changing 
their shapes; and the soli ton with large amplitude moves faster than 
the one with smaller amplitude. Thus, after a certain period of 
time, all these solitons interact and overlap in such a way that the 
initial waveform reoccurs. This remarkable phenomenon was first 
discovered in 1965 in the difference solution discussed above. This 
phenomenon has also been confirmed experimentally and found in 
several physical nonlinear systems. Because of this remarkable 
property, great efforts have been spent to devise analytical methods 
to predict the performance of these nonli ri 0 dr systGHis • A11 thGSG 

fascinating developments in the last ten years are summarized in a 
recent article by Scott[15], Although we have not considered any 
non-linear problem here, it can be said that the finite difference 
methods are a powerful tool in the analysis of linear and nonlinear 
dispersive waves. 



NORMALIZED DISTANCE 


Fig. 31 --Sol utions of the Korteweg-deV-iries eguation, 
5=0.022, at three different times. (tg=l/ii) 


The main thrust of this dissertation is in the development of 
some finite difference methods for linear dispersive waves. To this 
end we haveobtained simple difference equations for the solution of 
transient signal propagation in several kinds of stratified dispersive 
media. Wo demonstrated the numerical stability and accuracy of the 
difference equations, and described a number of practical techniques 
to greatly simplify the calculation procedure. The speed and ef- 
iciency with which time domain solutions may now be obtained using 
these methods lead us to expect many significant applications in 
the future. 




REFERENCES 


1. Sommerfel d , A., "About the Propagation of Light in Dispersive 
Media," Ann. Physik, Vol . 44, pp. 177-202, 1914. 

2. Brillouin, L., Wave Propagation and Group Velocity . New York: 
Academic Press, 1960. 

3. Wait, J. R. (Ed.), Electromagnetic Probing in Geophysics . 

Boulder, Colorado: The Golem Press, 1971. 

4. Gray, K. G. and Bowhill, S. A., "The Impulse Response of a Cold 
Stratified Plasma in the Presence of Collisions and a Vertical 
Magnetic Field by a Multiple Scattering Technique," Radio 
Science, Vol. 9, pp. 559-566, May 1974. 

5. Cronson, H, M. and Ross, 6. 6., "Current Status of Time-Domain 
Meterology in Material and Distributed Network Research," 

IEEE Trans. Instrum. Meas., Vol. IM-21, pp. 495-500, 

November 1972. 

6. Bowhill, S. A., "The Transient Response of the Ionosphere at 
Low Frequencies," The Physics of the Ionosphere, pp. 308-319, 
The Physical Society, London, 1954. 

7. Field, J. C., "Transient Propagation of Electromagnetic Waves 
in a Stratified Plasma," Radio Science, Vol. 6, pp. 503-510, 
April 1971. 

8. Courant, R., Friedrichs, K. 0., and Lewy, H., "On Partial 
Difference Equations of Mathematical Physics," Math. Ann. Vol. 
100, pp. 32-74, 1928. (English translation of original work, 
IBM J. Res. Devlop., Vol. 11, pp. 215-234, 1967). 

9. Milne, W. E., Numerical Solution of Differential Equations . 

John Wiley, New York, 1953, pp. 124-125. 

10. Chiu, S. S., "Difference Methods for Multiple Reflection of 
Elastic Stress Waves," J. Computational Phys., Vol. 6, 

pp. 17-28, January 1970. 

11. Casey, K. F., "Application of Hill's Functions to Problems of 
Propagation in Stratified Media," IEEE Trans. Antennas and 
Propagation, Vol. AP-20, pp. 368-374, May 1972. 


85 



86 


12. Frohlich, J., Theory of Dielectrics , Oxford University Press, 

1958. 

13. Zabusky, N, J. and Kruskal , M. D., "Interaction of Solitons in a 
Collisionless Plasma and the Recurrence of Initial States," 

Phys. Review Letters, Vol. 15, pp. 240-243, August 1965. 

14. Leibovich, S. and Seebass, R. A., Nonlinear Waves . Cornell 
University Press, 1974. 

15. Scott, A. C., Chu, F. Y. F, and McLaughlin, D. W. , "The Soliton: 

A New Concept in Applied Science," Proc. IEEE, Vol. 61, pp. 1443- 

16. Richtmyer, R. D. and Morton, K. W. , Difference Methods for 
Initial Value Problems , Second Edition, John Wiley (InterScience) , 
New York, 1967. 

17. Hildebrand, Finite Difference Methods and Simulations , New Jersey: 
Prentice-Hall, 1968. 

18. Forsthe, G. E. and Wasow, W. R., Finite Difference Methods for 
Partial Differential Equations . New York: John Wiley, 1960. 

19* Birkhoff, G. and Lynch, R. E., "Numerical Solution of the 
Telegraph and Related Equations," in Numerical Solution of 
Partial Differential Equations, Edited by J, H. Bramble, 

New York: Academic Press, 1966, pp. 289-315. 

20. Morse, P. and Feshbach, H., Methods of Theoretical Physics , 

New York: McGraw-Hill, 1953, pp. 139. 

21. Sommerfeld, A., Partial Differentia l Equations in Physics, 

New York: Academic Press, 1964, pp. 36-42. 

22. Stratton, J. A., Electromagnetic Theory, New York; McGraw-Hill. 

1941. 

23. Brekhovskikh, L. M., Waves in Layered Media . New York: 

Academic Press, 1960. 

24. Wait, J. R., Electromagnetic Waves in Stratified Media . Oxford, 
England: Pergamon Press, 1970. 

25. Wait, J. R., "Ground Wave Pulses," in Electromagnetic Probing 
in Geophysics, Edited by J. R. Wait, Boulder, Colorado: 

The Golem Press, 1971, pp. 349-360. 



87 


26. Sivaprasad, K. and Stotz, K. C., "Reflection of Electromagnetic 
Pulses from a Multilayered Medium," IEEE Trans. Geoscience 
Electronics, Vol . GE-11, pp. 161-164, July 1973. 

27. King, R, W. P. and Harrison, C. W, , Jr., "The Transmission of 
Electromagnetic Waves and Pulses into Earth," J. Applied Physics, 
Vol. 39, pp. 4444-4452, August 1968. 

28. Fuller, J. A. and Wait, J. R., "Electromagnetic Pulse Transmission 
in Homogeneous Dispersive Media," IEEE Trans. Antennas Propagat., 
Vol. AP-20, pp. 530-533, July 1972. 

29. Brillouin, L., "About the Propagation of Light in Dispersive 
Media," Ann. Physik, Vol. 44, pp. 203-240, 1914. 

30. Sommerfeld, A.. Optics . New York: Academic Press, 1954, pp. 88-91. 

31. Jordan, E. C. and Balmain, K. 6., Electromagnetic Waves and 
Radiating Systems . New Jersey: Prentice-Hall, 1968, pp. 729. 

32. Wait, J. R., "Propagation of Pulses in Dispersive Media," 

Radio Science, J. Research NBS/USNC-URSI, Vol. 690, pp. 1387- 
1401, November 1965. 

33. Haskell, R. E. and Case, C. T. , "Transient Signal Propagation 
in Lossless Isotropic Plasma," IEEE Trans. Antennas Propagat., 

Vol. AP-15, pp. 458-464, 1967. 

34. Gray, K. G., "An Exact Solution for the Impulse Response of a 
Uniform Plasma Half Space," IEEE Trans. Antennas Propagation, 

Vol. AP-22, No, 6, pp. 819-821, November 1974. 

35. McIntosh, R. E. and El-Khamy, S. E., "Compression of Transmitted 
Pulses in Plasmas," IEEE Trans. Antenas Propagation, Vol. AP-18, 
pp. 236-241, March 1970. 

36. Raptis, A. C., Mayhan, J. T. and Chen, C. S., "Pulse Compression 
in Bounded Dispersive Media," URSI/G-AP Symposium, Georgia 
Institute of Technology, Atlanta, Georgia, 10-13 June 1974. 

37. Wait, J. R., "Reflection of a Plane Transient Electromagnetic 
Wave from a Cold Lossless Plasma Slab," Radio Science, Vol. 4, 
pp. 401-405, April 1969. 

38. Hill, D. A. and Wait, J. R., "Reflection of Pulses from a 
Linearly Varying Ionosphere Model with a Vertical Magnetic 
Field," Radio Science, Vol. 6, pp. 933-937, November 1971. 



88 


39. Fante, R. L. and Taylor, R. L., "Effects of Losses on Transient 
Propagation in Dispersive Media - A System Study," IEEE Trans. 
Antennas Propagat., Vol . AP-21 , pp. 918, November 1973. (Also 
reports to which it refers) 

40. Fellner-Feldegg, H, , "The Measurement of Dielectrics in the Time 
Domain," J. Phys. Chem. , Vol. 73, pp. 616-623, March 1969. 

41. Fellner-Feldegg, H., "Permeability, Permittivity and Con- 
ductivity Measurements with Time Domain Reflectometry ," 
Hewlett-Packard Application Note 153, May 1972. 

42. Nicolson, A. M., Bennet, C. L., Lamensdorf, D. and Susman, L., 
"Application of Time-Domain Heterology to the Automation of 
Broad-Band Microwave Measurements," IEEE Trans. Microwave 
Theory Tech., Vol. MTT-20, pp. 3-9, January 1974. 

43. Fellner-Feldegg, H. and Barnett, E. F. , "Reflection of a Voltage 
Step from a Section of Transmission Line filled with a Polar 
Dielectric," J. Phys. Chem., Vol. 74, pp. 1962-1965, April 1970. 

44. Van Gemert, M. J. C. and De Graan, J. G., "Calculations on the 
Time Domain Reflection Behavior of a Transmission Line Partly 
Filled with a Polar Dielectric," Appl . Sci . Res., Vol. 26, 

pp. 1-17, May 1972. 

45. Antonucci , J. D., "An Artificial Transmission Line for Studies 
of Transient Propagation in Plasma Media," Air Force Cambridge 
Res. Labs., Bedford, Mass., AFCRL-72-0055, January 1972. 

46. Miura, M. R. , "The Korteweg-deVries Equation: A Model Equation 

for Nonlinear Waves," in Nonlinear Waves, edited by S. Leibovich 
and A. R. Seebass, pp. 212-234, Cornell University Press, 

Ithaca, New York, 1974. 

47. Hildebrand, F. B., Advanced Calculus for Applications , New Jersey; 
Prentice-Hall, 1965, pp. 27-28. 

48. Ito, M. , "Dispersion of a Step-Modulated Carrier Wave in Wave- 
guide," Proc. IEEE, Vol. 52, pp. 1250, October 1964. 

49. Proud, J. M., Tamarkin, P., and Kornhauser, E. T. , "Propagation 
of Sound Pulses in a Dispersive Medium," J. Acoustical Soc. 

Am., Vol. 28, No, 1, pp. 80-85, January 1956. 



APPENDIX A 
COMPUTER PROGRAMS 


Six computer programs which produce most of the results in 

this dissertation are collected here. 

1. Computer Program 1 calculates the reflected waveform of a 
unit-step sinusoidal wave normally incident on a specific 
inhomogeneous dielectric slab in [11]. Input data include: 

EPSR = 0.5 or 2.0; in Casey's paper (e^.=e.,/eT) 

NSTEP= number of time steps to be calculated'^ 

XL = normalized thickness L of the slab 
D = 1.0 if a transition layer, 

2.0 if a dielectric duct. 

2. Computer Program 2 calculates the reflection of a sine-squared 
pulse from a three-layered lossy dielectric medium[26]. Input 
data are: 

Z = thickness of the second layer 

NSTEP = number of time steps to be calculated 

DT = time step size (e.g.. At = 2.5 x 10-9 sec.) 

3. Computer Program 3 generates the results in Chapter IV - 
Sommerfeld and Brillouin's transient problem. Input data are 

NSTEP = number of time steps to be calculated 

WPDT = Up At 

WOWP = up/up 

PWO = p/up 

I0BSER= Zi/upAt + 1 

I0B2 = Z2/upAt + 1 

WWP = u/up. 

(Z 2 < Z-| where Z^ , Z 2 are spatial observation points) 

4- Computer Program 4 calculates the propagation and reflection of 
a transient unit-step, step carrier, and a chirp pulse in a 
homogeneous lossy plasma half-space or slab geometry. The 
input data are 

ICASE = 1 for propagation, 2 for reflection 
GWP — v/ Up 

NSTEP = number of time steps to be calculated 
WPDT = UpAt 
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Z = ojpz/c. In the case of propagation (ICASE=1) it 
is the observation point. In the case of re- 
flection (ICASE=2) it represents the thickness 
of the plasma slab. 

WWP * w/top, frequency of the step carrier if applicable 
IP =1 unit-step 

2 step carrier 

3 chirp pulse. 


5. Computer Program 5 computes the reflected waveform of a short 
gaussian pulse, a unit-step and a step carrier normally in- 
cident on an inhomogeneous lossy plasma (specifically the linear 
electron density profile[38]) or a homogeneous cold plasma 
half-space in which a longitudinal dc magnetic field is 
present* The input data are 


WPDTl = ojpAt if a homogeneous plasma, or a normalized 
factor for the incident short gaussian pulse. 
WCDT = oicAt 
GDT = vAt 
WWP — cij/ (jop 

NSTEP = number of time steps to be calculated 
IP =1 unit-step 

2 step carrier 

3 gaussian short pulse 

DT = 0 if a homogeneous half -space, or the time step 
size At otherwise. 


6. Computer Program 6 computes the reflected waveform of a 

unit-step signal normally incident to an ice layer on water 
in which water assumes a Debye model. The input data are 

NSTEPl = number of time steps to be calculated 
WODT = (DpAt (to=1/«o is the relaxation time) 

EPSl = dielectric constant of the second layer 
(ice in this case) 

EPS20 = ££ low frequency dielectric constant 
EPS2N = £„, high frequency dielectric constant 
CONDO = normalized low frequency conductivity of the 
Debye dielectric 
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COMPUTER PROGRAM 1 

REFLECTED UAUE FROM AN INHOMOGENEOUS DIELECTRIC SLAB 
REFERENCE: CASEY/ IEEE G-AP/ KAY 1972/ PP. 368-374 
DIMENSION F0C500)/ FI C 5 00 ) / F2 < 5 00 ) / EPS C 500 ) / DZ C5 00 ) / 

I VAFORM<500) 

DATA PI/PI2/Z0/FOLD/DZC 1) /3 . 1 4 1 59 / 6. 28 32/ 0 . 0/ 0 . 0/ 0 . 0/ 

READ<8/“) EPSR/ NSTEP/ XL/ D 

VRITE<6/9) NSTEP/ EPSR/ XL/ D 

FORMAK* NSTEP=* /I4// * EPSR =*/F4^2//* L' 

1' DCTRANSITION LAYER/D=U0; DUCT/ 0=2 • 0 > = • / F3 . I //) 
PROFILE OF THE INHOMOGENEOUS DIELECTRIC 
FDT=1 .0/36^0 
AA=4.0#FDT/XL 
DEL=(EPSH-1 •0)/(£PSR+l *0) 

EPS< U=1 -O-DEL 
DO 10 1=2/ 500 
EPS0 = 1 •0-DEL’!'COSCPI+Z0) 

DELZD=AA/SQRTCEFS0> 

EPSl =l .0-DEL*CCSCPI^CZ0+DELZD) ) 

EPSCI)=B.5+<EP50+EPS I ) 

DZ C I) =AA/S OPT C EPS C I ) ) 

Z0=Z0+DZCI ) 

1F(Z0.GT*D) GO TO 1 1 

IENDI=I 

IEND=IENDl+I 

VRITEC6/ 16> CSPSCI )/ 1=1/ lENDl) 

DO 12 1=1/ lENDl 
EPS< I >=SQKT(:EFS< I) ) 

CU SOURCE AT I=i; AIPVDIELECTHI C INTERFACE AT 1=2 

DO 13 X=l/IEND 

FOCI) =0-0 

FI C I ) =0.0 

FI C 1 )=SIN<PI2*FDT) 

FI (2)=2.0>t'EPS( I )>^F0( 1 )/CEPS( 1 ) +EPSC2) > 

VAFORMC 1 )=F1 C2) - F0<1> 

DO 15 N=2/ NSTEP 
F2< 1 )=SIN<N*PI2>«'FDT) 

IFCN-LT-IENDl ) M0VING=N+1 
DO 14 1=2/ MOVING 
T0TALF=F1C I-l ) 

IFCI •EO-2) T0TALF=F1 < 1 ) + C F0C 2 > - FOLD) 

T12=2,0/(EPSCI-1 )+EPS( I ) ) 

F2C I > = -F0C I >+Tl2*CEPSCI-l )^TOTALF+EPS( I)’«'F1 C I + I ) > 
F2CIEND)=F1(IENDI ) 

VAF0RMCN)=F2<2)-Fl( I ) 

FOLD=F0( I ) 

DO 15 1=1/ I END 
FBCI)=Fl(I) ‘ 

FI (I )=F2< 1 ) 

FORMATC* DIELECTRIC PROFI LE * // I 00 ( ! 0F6 • 3/ ) / ) 

URITEC6/ I7> CVAFORMCN) /N=l /NSTEP) 

FCP-MATC* CU REFLECTED WAVEFORM '// 1 00 Cl 0F6 . 3/ ) ) 

END 


PA6®J5 
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33 
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*i^Hi**m**ik** COMPUTE? P?OG?AM 2 ********** 

DIRECT TIME DOMAIN SOLUTION OF EMP REFLECTED FROM 
A THPEE-LAYERED LOSSY DIELECTRIC MEDIUM. 

REFERENCE; S I VAPRASAD . & STOTE ^ IEEE GEOSC. ELECTP.ONI CS > 
PP. 161-164-P JULY 1973. 

DIMENSION F0C500), F1C500)^ F2<500)^ UAFORM(500) 

DATA EPS l>EPS2-rEPS3/l .0, 7.0, 81.0/ 

DATA CONDI, C0ND2, COND3/0.0, 0.001, 0.001/ 

DATA PI, TAU, EPS0/3. 141592, 0.1E-6, 8.854E-12/ 
READ(8,“) Z, NSTEP, DT 
IEND = 2 NSTEP/2 

IENDI=IEND-1 

DZ=DT*(3 .0E+8/SQRTCEPS2) ) 

ID=Z/DZ 

I2=2+ID 

Ol^CONDl*DT/<2.0*EPS0*EPSl > 

O2=COND2>kdT/(2.0’^EPS0*E?S2) 

Q3=COKD3*DT/C2.0*EPS0*EPS3) 

C3=S0RT<EPS2/EP51 > 

C0=0.5*(C3^< 1 .0+Q2} + < I .0 + G 1 > > 

Cl =0.5*<C3>K< 1 .0-G2) + < 1 .0-Gl >) 

D0=l .0+Q2 
Dl=l .0-02 

E3=SQRTCEPS3/EPS2) 

Ee^^O.S’HCES^C 1 .3+C3)+C I .0+02) ) 

E1“0.5^<E3*< I .0-O3)+C 1 .0-02)) 

G0=( 1 .0+03) 

Gl“< 1 .0-03) 

DO 10 1=2, lEND 
F0(I >==0.0 
FI U )=0.0 

A SINE-SOUARE PULSE IS GENERATED AT 1=1 
F0C 1 )=0.0 

FI ( 1 )=(SIN(DT*PI/TAU) )**2 
FK2)=F0( D/C0 
VAFOPMC 1 )=FI (2)-F0C 1 ) 

EOLD=0.0 

DO 100 N=2, NSTEP 
T=N*DT 

IFCT .LE. TAU) F2 C 1 ) = ( S I N< T^PI /TAU ) ) ^+2 
IFCT .GT. TAU) F2C1>=0.0 
TCTALF=F1CI) + (F0(2)-EOLD) 

1=2 IS THE AIP/DIELECTPIC INTERFACE 
F2C2) = <-C1*F0(2>+C3»»^F1(:3)+7OTALF)/C0 
IFCN -LT. lENDl) MOVING=N+l 
IFCN .GT. lENDl) MOVI NG =MOUI NG - 1 
12 IS THE (2ND/3RD) LAYERS' INTERFACE 
DO 20 1=3, MOVING 

IFCI .LT. 12) F2CI )=< -Cl *F0CI ) +FI C I -I ) +F1 C I +1 ) >/D0 
IFC I .EO . 12) F2CI ) = < -El *FCCI ) + FK 2 - 1 ) +E3^F1 C I + 1 ) )/E0 
IFCI .GT. 12) F2C I ) =C -G 1 ‘♦‘F^C I ) +Fl ( I - 1 ) +F1 C I + 1 ) >/G0 
CONTINUE 
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56 

57 

58 

59 

60 
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62 
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64 


PICK UP THE REFLECTED UAVEFORM & STORE PREVIOUS FIELDS 
VAF0RMCN>=F2(2)“F1 C I > 

EOLD = F0( n 

DO 30 MOVING 

F0(I >=F1(I) 

30 Fia>=F2<I) 

100 CONTINUE 

Z=ID*DZ 

UHITE<6>2 1 ) NSTEP^Z^DT, CUAFOHM< N> , N= ! > NSTEP ) 

31 FORKATC* NSTEP=*>I4/^ ' Z DT =’>E10.3//> 

I’ SINE-SQUARED PULSE REFLECTED WAVEFORM • //^ 50 C 1 0F6 . 3/ ) > 

END 


1 NSTEP= 150 

2 Z * 4*9891 

3 DT = .220E -8 

4 

5 SINE-SQUARED PULSE REFLECTED WAVEFORM 

6 


7 

0m000 *•002 

-.009 

-.019 

8 

-•193 -*227 

-.262 

-.296 

9 

-.490 -.505 

-.5 16 

-.523 

10 

-.450 -.426 

-.399 

-.371 

11 

-.173 -.153 

-.138 

-.128 

12 

-.176 -.188 

-.199 

-.21 I 

13 

-.262 -.262 

-.261 

-.258 

U 

-.186 -.171 

-.155 

-.138 

15 

-.035 -.025 

-.016 

-.0 09 

16 

.008 

.0 10 

.012 

.0 14 

17 

.024 

.024 

.024 

.025 

18 

.016 

.017 

.015 

.013 

19 

.001 — .000 

-.001 

-.002 

20 

-.004 -.005 

-.005 

-.005 

21 

-.007 -.007 

-.007 

-.007 



.034 

-.053 

-.076 

-.101 

-.130 

-.161 

.330 

-.363 

-.394 

-.422 

-.446 

-.471 

525 

-.522 

-.516 

-.505 

-.490 

-.472 

341 

-.311 

-.261 

-.25 1 

-.223 

-.197 

125 

-.127 

-.134 

143 

-.154 

-.165 

222 

-.232 

-.241 

-.248 

-.255 

-.259 

253 

-.246 

-.237 

-.227 

-.214 

-.201 

122 

-.105 

-.089 

-.074 

-.060 

-.047 

004 

-.001 

.001 

.002 

.004 

.006 

016 

.017 

.019 

.020 

.022 

.023 

024 

.024 

.023 

.022 

.22 1 

.020 

01 1 

.009 

.008 

.00 6 

• 004 

.003 

003 

-.003 

-.003 

- .003 

- .004 

-.204 

005 

-.006 

-.006 

-.006 

-.006 

-.007 

007 

-.007 

-.007 

-.007 

-.007 

- .007 
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COMPUTE® PROGP.AM 3 *4t4«*****>it^ 

A NUMERICAL SOLUTION TO SOMMEPFELD'S ABOUT THE 
PROPAGATION OF LIGHT IN DISPERSIVE MEDIA’* — D* H. LAM 
DIMENSION F0C100 1 ) ^ FI ( 1 00 1 F2< KG 1 ) ^CC 1001 > ^SC 100 n ^ 
J UAFORMC 1000>^UAVE2e< 1G00> 

DATA II I^SS/UAF0?M< 1>jUAVE2SC 1 > /2 ^ 0 ♦ 0> 0 • 0^ 0 . 0/ 

READC8 / NSTEP^VPDT^'v 0 Iv^P>WUPj PUG > iOBSER^ 1052 

UR I TEC 6^ 10) N5T5P^VPDTiW0UP>UWP^PU0 

10 FORMAT C* NSTEP^* >14/^ * UFDT =*>FA.2/j* U0UP 

IF4^2/^* UUP PUG =*^F4*2//) 

I£ND=NSTEP/2+I0BS£R 
tend 1 = IEND-1 
NN1=NSTE?-I052 

IFCU0UP.EG*0.0.AND*PU0*EG*0.0> GO TO 11 

U0DT=V0VP*UPDT 

PDT=PU0*U0DT 

ADT=5QHTCV0DT^U0DT-?DT*PDT) 

A=UPDT*UPDT*2 # 0*PDT 

B=UPDT*UPDT* C ADT* ADT -PDT^PDT ) /ADT 

G=£XP(-?DT) 

CN=COS(ADT) 

SN=SINCADT) 

GN=G*CN 

11 Q2 = l •0+0^25=^W?DT>#'UFDT 
Q3=0*5^UPDT*U?DT 

DO 12 I-Ij. I end 
CC n=^0*G 

S< I 

F0CI )=0K 
FlCI)=0*0 

12 F2CI)=0.0 
UDTsUUpsKV’PDT 
FI ( I >=SINCUDT> 

F1<2)=F0< I >/Q2 
DO 21 N=2^ NSTEP 
F2C 1 )=SINCN*UDT) 

IFCN -LT* I END I ) MCVING=W+1 
IFCN •GT. lENDi) MO VING^^MOVI NG - 1 
IFCN .GT. NNl > Ill-IIl+l 
DO 20 1=111, MOVING 

IFCU0U?.EOKK.AND»PU0*EO*0*0> GO TO 2C 
SOLD=SCI ) 

SC I )=G’»'CCN>S C I ) +5N*C( I ) +0 • 5* SN*F0 C I) ) 

CC 1 >=G*CCN*CC I ) -SN*SOLD)+0*5»ic( n C I ) +GN*FGC I ) ) 

SS=A*CC I )+3»t«SC I > 

20 F2CI)=-Fecn+CFl CI-1 )+FlCI^l)-C3><‘Fl C I) + SS)/Q2 
U A FO P.M C N ) = F2 ( I C 3 5 EP ) 

UAVE2G CN') =F2C I CF2) 

DC 21 1=1, MOVING 
F0CI)=F1C1) 

21 F1CI)=F2CI) 

URITEC6,22) I OES EP, CUAFCRMC N) , N= 1 , NSTEP ) 

22 FOPMATC* UAFORM SEEN AT 1 0E5EP= ’ , I A/, I 00 C 1 0F6 • 3/ ) ) 
END 
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21 


22 


23 


2A 

S 

25 


26 


27 


23 

33 

29 


30 


31 


32 


33 


3A 


35 


36 

1 1 

37 
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30 

SI 



CO:<PUTEF. PPOHRAM A + 

THIS PP.0G?.A>1 calculates THE PHOPAGATIOM AMD PEFLECTIOM 
OF A LIMIT- STEP# CU, PHASE- VAH I AT I OM PULSE IM AM 
ISOTROPIC PLAS-MA HALF- SPACE OH SLAB GEOMETP.f. 

DI MEMSIOM FO< 800)^ r H 300)^ F?5 C 3 00 ) ^ SU.M ( 3 00 ) j UA FORM C 800) 
COH^MOM UDT^ WPDT> IP 

PEADC8#-) ICASE# GW?, MSTEPj WPDT* 6, WWP, IP 

WD T = W v; P * UP D T 

lD=<(i+O.OOOOU/UPDT 

IS=3+1D 

IQHSEP=^/WPDT + I 

I F( I CASE .EO* a> I03SER=2 

iemd=iobse;h + CMSTEP- I0BSER>/3 1 

IEMDUIEMD-1 

GDT^GUFtUPDT 

EaP1=EHPC-GDT) 

3A= GDT^f^UPDT* >7PD T 

QA5=0*5+0A 

Q3=0.5*UPDr4t:OPDT 

035=03/2*0 

Ql= I *0+035/2*0 

02=1.0+03/2*0 

VRITE(6jH> gup* UPDT^ lEMD^ MSTEP> i 

FOPMATC GDT/VPDT = *,F7.3/>* UPDT = ’^F7*3/> 

2’ « GRID POIMT USED = 9 TI.<E STEPS HUM = SIA//> 

3* UP-^ii/C = ’ j F5. 1//) 

UP I TE C 8 3 > 1 C A SEj i P > UUP > 1 0 BSEP 

FOrtHAK* ICAS£(PS0PAGAT10M= i; HEFLECT I 0 M= 2 ) = Sip/* 

1’ IPCUMI T-SrEP= 1 ; CU UAVE=2> = SIP/> 

2’ l\^L^(FHE0lj?:MCf OF CU SOURCE) = %FS*2//, 

3* UAFORM OBSERVED AT GRID POIMT = %I3//> 

DO U 1=2# lEMD 
SlIMC I > = 0*0 
F0Cn = 0.0 
FI Cl ) = 0*0 
F2( I )=0*0 
FOCI )= PULS ECO) 

FI C \ )= PULSE C 1 > 

GO TOC 12# 13)> I CASE 
FI C9)=F0C I )/Ql 
yAFO:PH( n = FI C 2 )-pulsec o) 

EIMC=0*0 
GO 10 1 A 
FI (2)=F0( 1 )/09 
UAFOHMC 1 ) = F1 C lOFSF":) 

DO 21 M=2# MSTEP 
F2C 1 ) = P :lsec M) 

IFCM *LT* lEMDl) M0VIMG=M+1 
GO TOC 30# 31)# I CASE 
I FC M* GT. I EMDl ) H )VI MG=HOVl MG- 1 
no P? I=P# /iOVIMU 


0? 
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I P< 0^1 •GT« ) SU.^i ( I ) I ) ♦SXP 1 + 5-^'< F.X9 I ^ F-:^ ( I > + F U I > ) 

53 21 FP< n =-FC^( I ) + ( F U I - 1) + FU I + I ) i )+Q^i+SU/*C \ }^ /Q2 

54 vJAFCH-M(iNi) =F^?C I C3S£R) 

55 GO TO 32 

56 C I=5’:> IS THE AlH/PLASiOA I.MTFRFACE 

57 31 S0URCF=F1 ( 1) + ( F0( P > - Ki :'JC ) 

55 i=2 

59 1 F< Q4 • G T • • 0 ) SUM C 1 > = S U/i C I ) * FXP 1 + 'I • 5^ < FXP I + F3 ( I) + F U 1 > > 

60 FP( I >=-F0< I SOUHCE+Fl ( I + I ) ^ 33 5* FI C I ) + 04 5+SUM( 1 ) > / Q 1 

61 DO 26 I=3> MOyiiVG 

6S 1F(Q4^GT.0* *A*VD* I*LF.IP> SUM C 1 ) =S UM ( I ) + RXP 1 ^ 

63 1 0.5*( KXP1 + F0C 1 ) + r IC I )) 

64 IF( I -LT* 1 P) FP< I ) = -F3C n + ( F) C 1 - I ) + F U 1 + n -33*F1 ( 1 > 

65 I +Q4i‘SaMC I ) ) /3P 

66 I F( I .KQ- I P) FPC n =-F0( I ) + ( FI ( 1 - I > + FI ( 1 + n -Q35^ FlC r ) 

67 1 •►Q45<‘SJM( 1 ) ) /Q1 

68 26 IF(I*GT*1P) Fp ( I ) = F 1 C I - 1 ) 

69 WAFOrtM<v>J ) =F^?( P) -Ft C 1 ) 

70 EIMC=F0tl) 

71 32 CO.\)riMUF 

72 C SIMULATIOlM 5/14/74 


7 3 F2il Fi^JD ) =F 1 ( I FiM D- I > /O 1 

74 F0( IFi\ID)-Fl ( IR.viD) 

75 Fl(lF.i\in)=FPCIFMD) 

76 DO PI 1 = 1# MOO I MG 

77 FPKI>-Pl<n 

iH ai Fi(n = Fp<j) 

7 9 *^rtITFC6^pP) ( ‘a’AFORMCM) >M=:1 j MSTFP> 

HPS 22 FORMATC i:i :u nF6* 3/) > 

H 1 FMD 

H2 FUMCTIOM PULSF(M) 

83 COMMOM '.^D r> JPDT* Ip 

HA IFtlP 3) GO TO 

85 IFCIP .FQ* 1) ?rjLSF=1.0 

86 IK(IP .FQ. P) PULSF-SKNJUvUr*M> 

87 RETURN 

88 100 IFCM.nr.^) GO TO 101 

89 P=^wP 

90 H=4*n 

9 1 Q= f/SQPTC I •0-pi'?> 

92 S1^.-«PDT/:;) 

93 rOH-0- {t :i/SQHT( ri=»‘r^-p4<P) 

94 IFl^PDT •MR* Ma;«t= 1 + TCR/a/PDT 

95 ini PJLSF=0.0 

96 IFCM .GT* .'jOP) HFTOPM 

97 SM-M^Sl 

9 8 CFi‘A = iJ* ( P“ SGRTC ♦ SM- S->< ♦ P * P > ^ 

99 PULSF-5 [t\/< nPTA ) 

103 r.FIUPM 

10 1 F.>^n 
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;k4c*>k**>K**A COMPUTER PROGRAM 5 

■THIS PROGRAM COMPUTES THE RETLECTED WAVEFORM OF A SIGNAL 
FROM A HOMOGENEOUS AMD I NHOMOGEMEOUS PLASMA HALF-SPACE 
WHERE A VERTICAL MAGNETIC FIELD IS ASSUMED 
. DIMENSION EX0C510)>SXK510)>EX2C510);r5X(510>^SY(51O)> 

1 EY0(5 105^EYI C510)^EY2<51{3)>CXC5 10) / CYC510), 

2 REFEXC I00Z)^REFEYC I 020 ) > VPDT < 5 1 0 ) 

COMMON IP^ VDT^ VPDT 1 

, DATA EX0jEXl>EY0,EYl^CX^CY^SX-r5Y>EINC/510*0*0^510’i‘0.0^ 

I 5 1 5 10*0 • 0^ 5 1 0*0 *0^ 5 1 0*^^0j 510*0*0>510^0*0>0,0/ 
READ(8^-> UPDT!^ WCDT^ GDT^ WI7F/ NSTEP^ IP, DT 
WR I TEC 6^4) DT^WPDT 1 ^ WCDT> G DT^ UUP^NSTE?^ I ? 

FORMATC’ DT = '^F6*3,» NANOSECONDS’/-.^ WPDT 1 = ’ F6« 4// 

1’ UCDT=’^F6*4/> * GDT=’ ^F6-4/-. * TOP=’>F6.4/> 

2’ NSTEP=* , 14/^ * IPCUNIT-STEP=I-rCW WAVE = 2-. IMPULSE =^3 ) = ’ 
3-rI2//,.* REFLECTED WAVEFORMS; EX 6 , SYV) 

FOR IMPULSE SIGNAL : WPDTl-0.02 WITH 600 TIME STEPS 

IEND=2'fNSTEP/2 

IEND1=^I£HD-1 

DZ=0.3*DT 

LINEAR ELECTRON DENSITY PROFILE 
DO 5 1=2^ lEND 
Z=<I-l .5>*DZ 

WPDTC X>=0*000 14142*DT*SQRT<2 > 

WPDTC O«0.0 
IF(DT*NE.0*0) GO TO 7 
DO 6 1=1, LEND 
WPDTCI )=WPDT1 
AA=^WPDT< 1 ) WPDT C I ) 

DO 8 I ==2^ lEND 
BBssWPDTC I ) WPDTC I ) 

WPDTC I >=0.25 >««CAA+EB) 

AA=BB 

C=COSCWCDT) 

S=$INCWCDT) 

GG=EXPC-GDT) 

GGS=GG=^S 

GGC=GG>!'C 

DDI =GG*CWCDT*$+GDT*C) 

DD2=GG^CWCDT=<«C-GDT*£) 

WDT=WPWPDT1 
EX0C 1 )=PULSE(0> 

SXl C ! )=?ULSEC I > 

EXl(2)=EX3CI )/C I .G+0.5WPDT<2>> 

REFEXC 1 )=EX1 C2)-EX0( 1) 

REFEYC 1 >=EY1(2> 
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FINITE DIFFERENCE SOLUTION STARTS HERE 
DO 1 1 N=2/ NSTEP 
EX2C1 )=PULSE(N) 

EY2( i >^EYI C2> 

IFCN -LT. lENDl )M0VING=N+1 
IFCN *GT. lEWDl )MOVING=NOVING-l 
DO 10 1 ~ 2 j moving 
01 =yPDT( I) 

Q2 = I •0'f0.5*Qi 
U2=2p0>t^Ql 

SOLD=SX<n 

SXC I) = CGGC*SX( I )+GGS*CXC I > ) +0*5*CGDT*EX I ( I )+DDl*EX0CI ) ) 
CXd ) = CGGC*CXCn -GGS*SOLD)+0*5=i=CUCDT*EXl (I )+DD2*EX0( I > ) 
SOLD=SYU ) 

SYC I )-CGGC*SYC I > +GGS''J'CYC I ) > +0 • 5* C GDT + EY I C I ) +DD1 ’J‘EY0< I ) ) 
CYC I )==<GGC*CY( I ) -GGS + SOLD) +0 • 5*CUCDT*£Y1 ( I )+DD2*EY0C I ) ) 
S0RCEX=EX1 <1-1 ) 

IFCI -EQ.2) SOP.CEX^EXl ( 1 ) + C EX0( 2 ) -E INC ) 

SSX=SX(I)-CY(I) 

EX2C I >=-EX0< I) + (EX1 < I + D+SOPCEX -0 I *EX I < I ) +U2*S SX ) /G2 
SSY=CX( I )+SY(I ) 

EY2C I)=-EY0(X)+CEY1< I + l)+£Yl CI-l)-Gl*EYI (I ) +V2=^SSY) / 02 
CONTINUE 

PICK UP THE REFLECTED WAVEFORM <& UPDATE THE FIELDS 
REFEXCN>=EX2<2>-EX! C 1 ) 

REFEYCN)=EY2(2) 

EINC=EX0C 1 ) 

DO 11 1=1, MOVING 
EX0CI)=EX1 (I ) 

EXICI )=EX2(1 ) 

EY0CI )=EY1 C 1 ) 

EY1CI)=EY2(I) 

UR1TEC6, 12) <REFEXCN>,N=1, NSTEP) 

FORMAT<100ClGF7.4/)> 

URITE(6,12) (HEFEYCN),N=1 , NSTEP) 

END 

FUNCTION PULSECN) 

COMMON IP, VDT, UPDT I 

DATA K, N0,NN2, TP/4, U, 22, 6*2332/ 

GAUSSIAN PULSE PARAMETERS: N0=K^2*6; NN2=2*N0 

PULSE=0.0 

IFCIP -EQ. n PULSE=U0 
IFCIP *E0- 2) PULSE=SINCN’^'UDT> 

I F< I P.EQ.3 *AND*N.LE.NN2> 

I ?ULSE=EX?C - C CN-NS) /CK^ 1 *4 I 42 ) > ) /SCRTCTP* C K*V?DT 1 ) 

RETURN 
END 
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COMPUTER PROGRAM 6 ♦Hc****%**><5 
THIS PROGRAM COMPUTES THE REFLECTED UAVEFORM OF A 
UNIT^STEF SIGNAL NORMALLY INCIDENT TO A THREE -LA’^ETP£Q 
DISPERSIVE MEDIUM^ AI H/I CE/UATEP^ IN WHICH WATE^ 

ASSUMES A DEBYE MODEL • 

DIMENSION F0<510)-fFl C510)jF2C5I0>^ SUMC 510) ^ WAFORMC 2000 ) 
READCe^-) N5TEPI> W0 DTj EPSl^ EPS20^ EPS2N^ COND0 
IEND=2+NSTEPl/2 
1END1=IEND-I 

VRITE<6^ 102) NSTEPl^ W0DT^ EFS1> EPS20^ EPS2N^ COND0 
102 FORMATC ^ OF TIME STEPS RUN = ’>15/, • W0DT = ’,F3.5/, 
DIELECTRIC C0NSTANTC2MD LAYER) ^ ’>F5.2/, 

A’ LOW & HIGH FREQ DIELECTRIC CONSTS* OF WATE^=^ ’ , 2F9 . 3/, 
2* COND0(COND*TAU0/EPSAIR) OF WATER *,F10*2//) 
EXP1=EXPC-W0DT) 

EPS8=U0 

CONDR=COND0/EFS2N+ ( EPS 20 -EPS 2N) /EPS 2N 
01 = 1# 0+0* S^i'CONDR+WODT 
02=1*0-0* 5^CONDR*W0DT 
Q3=W0DT>i<W0DT+<EPS20-EFS2N>/EPS2W 
Q4 = Q3’«'W0DT 

EPS21=SQRT<EPS 1/EPS2N) 

Q5=CSPS2I+Q1 >/2.0 
Q6=(EFS21+O2)/2*0 
DO M 1=2, lEND 
SUMC I ) =0*0 

reel >= 0.0 

11 FIC1)=0.0 

REFLECTION & TRANSMISSION COEFF* AT THE AIR/ ICE BOUNDA®^ 
DUM=SQRTCS?S0>+SORTCEPS I ) 

T0 I =2 .0*SQRTCSPS0) /DUM 

TI0=2.0*SORT(EPS1 )/DUH 

R01=<SORTCEPS0>“SORTCSPS1 ) )/DUM 

INCIDENT FIELD TO THE WATER SURFACE 

F0( I ) = 1 *0*T6I 

Fl< I ) = 1 *0*TS1 

FI C2)=EPS2I*F0C 1)/Q5 

WAFORMC 1 )=R0i + CFl C2)-T01)=!'T10 

FINITE DIFFERENCE SOLUTIONS START HERE 

DO 200 N=2, NSTEPl 

F2C I )=FI C2) 

1-2 IS THE ICE/WATER INTERFACE 

1=2 

SUMC I ) =SUMCI >*EXP1 + <F0( 1 )>»<EX?1 +F1 C 1 ) )/2 -0 
F2 C I ) = C -06^FGC I) +F1 Cl + 1 ) +EPS2 I 1 C I-l)+G.5*(G3*FKI)- 
I 04+SUMC I > ) )/05 

EFFICIENT MOVING TIME WINDOW CALCULATION‘S ADO^=’T^D 
IFCM'.LT* lENDl) M0VING=N+1 
IFCN *GT. lENDl) MOVING =N0 VI NG - 1 
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DO 20 1=3^ MOVING 

5UM( I > =SUMC I )*EX?1 + C F0C n^EXPl + FI ( ! > )/2.0 

F2CI ) = < -'Q2’^<F0C I ) +Fi C I - 1 )+Fl ( I + 1 > +23*F1 C I ) “04l*5UM( I ) ) /G 1 

20 CONTINUE 

UPDATE THE FIELDS AT A NEU TIME STEPS 

DO 21 1=^1-. MOVING 

F0(I}:=FICI) 

21 FI(I/=F2(I) 

PICK U? THE REFLECTED UAFOPM AT THE AIR/ ICE INTERFACE 
VAFOm'KN)=r!01+CFl<2>-T0n^Tl0 
200 CONTINUE 

VP. I TE C 6.* 23) R0 ! 

23 FORMAT (’ REFLECTION OF THE AIR/ ICE BOUNDARY = *.»FIS.4//> 
WRITE <6/ 105) 

105 FORMATS ’ REFLECTED WAFORM OF A UNIT-STEP SIGNAL FROM’/ 

1' THE ICE/WATER INTERFACE SEEN IN FREE SPACE’//) 

WR I TE < 6 / 2 2 ) CW AFORM NSTEP 1 ) 

22 FORMATC^00( 10F8.4/) > 

END 



